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ABSTRACT 


Inversion  of  Laplace  transforms  has  been  accomplished 
by  a  numerical  integration  along  appropriate  paths  in  the 
complex  plane.  Two  general  procedures  have  been  used. 

The  simpler  and  more  economical  employs  a  simple  path,  such 
as  a  parabola,  which  bends  to  the  left.  Accuracy  is  main¬ 
tained  by  monitoring  the  oscillation  of  the  integrand.  A 
second  method  employs  a  steepest  descent  contour. 
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I .  INTRODUCTION 


If  F(t)  is  an  integrable  function  of  the  real  variable 
t,  the  function  f(s)  defined  by  the  integral 

oo 

f ( s )  =  j  e"st  F(t)dt  (1) 

o 

is  called  the  Laplace  transformation,  or  the  Laplace  trans¬ 
form,  of  F(t)  and  is  frequently  indicated  by  the  notation 

f(s)  »£[F(t)]  (2) 

The  Laplace  transformation  is  a  linear  integral  trans¬ 
formation  which  is  widely  used  in  applied  mathematics  and 
technology.  A  typical  application  is  as  follows.  An  un¬ 
known  function  F(t)  satisfies  a  certain  differential  equa¬ 
tion  and  specified  end  conditions.  Employing  theorems 
applicable  to  the  Laplace  transform,  a  function  f(s)  is 
determined  as  the  solution  of  an  algebraic  equation  which 
corresponds,  in  an  appropriate  fashion,  to  the  differential 
equation  and  end  conditions  of  F(t).  Then  the  problem  is 
reduced  to  the  following  task:  having  established  f(s), 
it  is  necessary  to  determine  the  function  F(t)  for  which 


This  problem  is  called  finding  the  inverse  (Laplace) 


transform  of  the  function  f(s).  This  is  customarily 
written  as 

F(t)  =<£"V(s)]  (3) 

It  may  be  shown  that  the  following  inversion  integral 


F(t)  =  J  ^  f(s)ds  (4) 

c 

accomplishes  the  desired  inversion.  In  this  expression 
s  is  the  complex  variable 

s  =  x  +  iy  (5) 

and  the  integration  is  along  a  suitable  path  in  the  complex 
s-plane. 

If  in  fact  there  is  a  function  F(t)  such  that  the  given 
f(s)  -*.[F(t  )],  then  the  function  determined  by  equation 
(4)  is,  indeed,  the  function  F(t).  The  result  is  unique 
in  the  following  sense.  If 


and 


f ( s )  =  J  e-st  F1(t)dt 

o 

F2(t>  "  35T  f  eSt  '<s)ds 

c 


(6) 


(7) 


then  F^  ( t )  and  F2(t)  differ  only  on  a  set  of  Lebesgue 
measure  zero. 

In  this  thesis  numerical  methods  will  be  utilized  to 
perform  the  integration  of  equation  (4)  and  functions  for 
which  the  inversion  is  successful  will  correspond  to 
functions  F(t)  which  are  at  least  piecewise  continuous, 
having  only  isolated  discontinuities,  if  any  are  present 
at  all.  If  F(t)  is  discontinuous  for  t  =  a,  then  the 
recovered  function  satisfies  the  condition 

F( a)  =  |  Lim[F(a+e )  +  F(a-e)]  (8) 

e  -  0 

The  path  of  integration,  indicated  by  C  in  equation 
(4)  is  frequently  described  as  an  infinitely  long,  vertical 
straight  line  x  =  constant,  to  the  right  of  any  singulari¬ 
ties  of  f(s).  (A  theorem  of  complex  variables  shows  that 
f(s)  is  analytic  in  its  region  of  convergence,  a  right-hand 
plane,  except  for  isolated  singularities.)  This  straight 
vertical  contour  is  frequently  called  the  Bromwich  contour. 

Usually  inversion  of  a  Laplace  transform,  that  is  find¬ 
ing  F(t)  =  [  f(s ) ]  ,  is  accomplished  by  use  of  a  table  of 

transform  pairs.  One  of  the  most  complete  of  these  is  that 
of  Roberts  and  Kaufman  [14].  Other  such  tables  are  listed 
as  references  [5],  [6],  [7],  and  [13], 

Several  theorems  concerning  the  Laplace  transform  extend 
the  usefulness  of  such  transformation  pair  tables.  For 
example 
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<£"1[f1(s)  +  f2(s)]  =od'1rf1(s)]  +aC’1[f2(s)] 

(9) 

However,  it  is  not  an  infrequent  occurrence  that  one  is 
called  upon  to  find  F(t)  in  a  case  where  the  tables  and 
the  theorems  are  of  no  direct  help.  There  has  developed 
a  considerable  literature  concerning  this  problem. 

One  obvious  method  is  to  attempt  to  expand  the  given 
f(s)  as  a  series  of  functions  fi(s)for  which  the  inverses 
are  known.  Then  using  the  linear  property  indicated  in 
equation  (9),  the  result  may  be  expressed  as 

<£_1[f(s)]  =6C"1[Za.f.(s)]  =  Za.«dr1[f.(s)]  (10) 

The  integration  indicated  in  equation  (4)  can  also  be 
performed  by  various  analytical  and  approximate  methods. 

In  particular,  one  can  perform  the  integration  along  the 
Bromwich  contour  by  a  numerical  procedure.  Much  of  the 
literature  is  devoted  to  such  methods. 

In  general  these  methods  amount  to  the  following.  A 
finite  number  of  points  along  the  Bromwich  contour  are 
determined  according  to  some  law  and  the  integrand  is 
evaluated  at  such  points.  There  are  associated  weight  fac¬ 
tors  and  the  integral  is  approximated  as  a  sum  of  products 
of  integrand  values  times  weight  factors.  One  of  the  most 
widely  used  of  such  methods  is  described  by  Salzer  [15], 
[16],  and  [17].  Other  procedures  of  the  type  are  treated 
in  [18],  [20],  [2],  [8],  [9],  [10],  and  [21]. 


It  is  interesting  to  note  Inat  the  same  classical  poly¬ 


nomials  which  are  useful  in  the  expansions  indicated  by 
equation  (10)  are  also  encountered  in  locating  points  for 
evaluation  using  a  numerical  integration  scheme  so  that 
the  final  algorithm  may  be  the  same  even  though  the  original 
motivation  was  quite  different. 

Other  methods  which  do  not  fall  clearly  into  the 
classes  of  methods  discussed  above  have  been  presented  by 
Widder  [22],  and  Bellman,  Kalaba,  and  Lockett  [3], 

These  approximate  methods  have  not  always  been  success¬ 
ful.  Hiep  [12]  employed  Salzer's  method  for  numerical  in¬ 
version  on  a  problem  of  heat  transfer  in  porous  media  and 
found  it  difficult  to  obtain  results  which  did  not  exhibit 
physically  incorrect  behavior,  e.g.,  in  the  vicinity  of  a 
point  of  engineering  interest,  the  temperature  of  the 
medium  fell  below  sink  temperature  and  rose  above  source 
temperature.  In  similar  work  on  a  conjugated  heat  transfer 
problem  Zargary  [23]  encountered  the  same  difficulty,  and, 
after  devoting  considerable  time  and  attention  to  the  prob¬ 
lem,  was  forced  to  abandon  the  use  of  numerical  inversion. 

To  overcome  such  difficulties  in  problems  of  these 
kinds,  we  have  reverted  to  direct  numerical  integration  of 
equation  (4)  along  appropriate  contours  in  the  s-plane. 

Our  procedures,  which  appear  thus  far  to  be  efficient  and 
reliable,  are  described  in  what  follows. 

U 


utmm 


II.  DISTORTED  CONTOURS 

We  have  found  that  great  advantage  accrues  from  distort¬ 
ing  the  Bromwich  contour  in  ways  which  will  subsequently 
be  described. 

According  to  theorems  concerning  analytic  functions  of 
a  complex  variable,  the  Bromwich  contour  can  be  distorted 
arbitrarily  as  long  as  the  process  of  distortion  does  not 
cause  the  contour  C  to  cross  over  a  singularity  of  f(s)  and 
as  long  as  the  two  ends  of  C  extend  to  infinity  in  a  manner 
which  preserves  the  convergence  of  the  integral  in  equation 
(4). 

Thus,  with  the  supposition  that  f(s)  is  analytic,  with 
the  exception  of  isolated  singularities,  and  is  real  valued 
for  real  s,  let  g(s)  represent  the  integrand 

g ( s )  »  eSt  f ( s )  *  u(s)  +  iv(s)  (11) 

where  u(s)  and  v(-e)  are  the  real  and  imaginary  harmonic 

components  of  g(s).  Furthermore,  assume  that  the  contour  C 

*  *  * 

consists  of  the  two  symmetrical  parts  C  and  C  as  shown 

in  Figure  1.  Consideration  of  the  contributions  to  ^"1[f(s)] 

*  ** 

over  arcs  ds-^  on  C  and  ds2  on  C  yields  the  following: 


dx2  ^-dx^ 

dy2  =  dyx 
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(12) 


FIGURE  1.  Distorted  Bromwich  Contour.  Contour  is 
distorted  into  C*  and  its  reflection  C** 


Thus, 

^[g1ds1  +  g2ds2]  =  |[v1dx1  +  UldYl]  (13) 

and,  therefore, 


F(t) 


[vdx  +  udy] 


(14) 


Equation  (14)  is,  of  course,  further  simplified  if  C* 
is  selected  to  be  a  path  along  which  v*0.  Additional 
theoretical  reasons  exist  for  such  a  selection.  Such  a  path 
will  subsequently  be  referred  to  as  a  path  of  steepest 
descent.  Numerical  integration  along  such  a  path,  or  a 
close  facsimile  thereof,  was  done  by  Esc h  [11]  in  1957  and 
Carrier,  Krook,  and  Pearson  [4]  recommended  further 


exploration  of  this  procedure.  However,  no  further  such  in¬ 
vestigation  seems  to  have  been  done  until  the  present  thesis. 

In  Chapter  V,  a  description  of  an  algorithm  for  performing 

* 

numerical  integration  along  a  path  C  ,  for  which  v=0  over 

a  substantial  portion  of  its  range,  will  be  presented. 

Almost  by  accidental  discovery,  the  present  work  has 

revealed  that  it  is  computationally  more  efficient  simply 

* 

to  make  an  arbitrary  choice  for  the  contour  C  which  satis¬ 
fies  the  following  criteria: 

* 

(1)  C  is  a  smooth  continuous  curve  of  the  form 

s(p)  -  x(p)  +  iy(p)  (15) 

where  p  is  a  real  parameter  0  £  p  £  00 . 

(2)  C*  lies  to  the  right  of  all  singularities  of  f(s). 

* 

(3)  C  does  not  extend  sufficiently  far  to  the  right 

to  cause  computational  difficulties  due  to  the 

factor  est. 

* 

(4)  C  approaches  infinity  in  such  a  way  that  x 
approaches  -®  and  y  approaches  ®. 

(5)  Oscillation  of  the  integrand  is  not  excessive. 
(Otherwise  there  may  be  loss  of  accuracy  due  to 
the  sampling  rate  of  the  integration  algorithm 
seriously  mismatching  variations  in  the  integrand 
and/or  positive  and  negative  contributions  nearly 
cancelling  each  other.) 

Such  a  contour  will  be  referred  to  as  a  simple  parameterized 
contour.  Discussion  of  how  to  assure  satisfying  criteria 


15 


(2),  (3),  and  (5)  will  be  deferred  until  after  the  develop¬ 
ment  of  the  actual  algorithm. 

The  integral  shown  in  equation  (14)  is  approximated  by 
the  simple  trapezoidal  sum 

M 

F(t)  ■  {[v(sk)+v(iWi(xk+i  -  xkJ 

k=l 

+  [u(sk)+u(sk+1)][yk+1-  yk]>  (16) 

where 

Sk  =  xk  +  iyk  =  x(pk)+  iy(pk)  (17) 

The  arbitrarily  chosen  functions  x(p)  and  y(p)  are  monotoni- 
cally  increasing  and  p^,  Pg,...,  Pn  is  an  increasing 
sequence. 

A  very  simple  example,  but  one  which  has  been  employed 
very  successfully,  is  given  by  the  following  equations 

x  *  A  -  Bp2 

y  -  P  (18) 

p  *  0  ,  A  ,  2A ,  .  .  . 

The  real  number  A  is  selected  so  as  to  satisfy  criteria  (2) 
and  (3)  above.  The  real  number  B  has,  usually,  been 
assigned  a  value  of  one. 

The  sum  includes  only  a  finite  number  of  terms.  It  is 
necessary  to  provide  an  appropriate  criterion  for  terminat¬ 
ing  the  process  of  summation.  This  will  be  discussed  later. 


In  order  to  monitor  the  possibility  of  inaccurate 
evaluations  resulting  from  oscillations,  in  addition  to  the 
sum  F(t)  given  by  equation  (16),  the  sum  G(t)  given  by 

M 

0<t)  -  sr  2  {[v(sk+i)-v(sk)i[xk+i-xkI 
k*l 

+  [•(Vi),Iisk)1[Vi-,k11  (19) 

is  also  obtained.  This  is  the  sum  of  the  absolute  values 
of  the  addends  to  the  sum  (16).  Additionally,  the  number 
of  times  successive  addends  are  of  opposite  signs  is 
recorded. 

The  present  numerical  experiments,  to  be  described 

later,  were  generally  so  successful  that  the  effect  of 

oscillation  could  not  be  seen.  However,  in  some  cases 

* 

where  poor  contours  C  were  deliberately  chosen,  it  was 
found  that  inaccurate  results  were  obtained  if  G(t)  was 

g 

many  orders  of  magnitude  (e.g.,  10  )  times  as  great  as  F(t) 


III.  THE  COMPUTER  PROGRAM  AND 
ITS  PRINCIPAL  SUBROUTINES 


Previous  discussion  proposed  two  alternatives  in  terms 
of  advantageous  distortion  of  the  Bromwich  contour  upon 
which  to  perform  the  numerical  integration  of  equation  (16). 
The  simple  parameterized  curve  has  been  experimentally  found 
to  be  the  more  efficient  of  the  two  methods  and  its  algo¬ 
rithm  will  be  developed  at  this  time.  The  steepest  descent 
procedure  will  be  treated  in  Chapter  V. 

The  computer  program  and  subprograms  which  implement 
the  simple  procedure  developed  in  the  preceding  chapter  are 
all  written  in  the  FORTRAN  language,  using  double  precision 
arithmetic.  They  have  been  tested  and  debugged  on  an  IBM 
360/67.  For  machines  having  a  larger  mantissa,  single  pre¬ 
cision  arithmetic  might  prove  to  be  satisfactory. 

SUBROUTINE  VALUE  is  a  user  supplied  routine  which  imple¬ 
ments  computation  of  the  g(s)  given  by  equation  (11)  for  the 
desired  transform  f(s).  This  subroutine  calculates  the  u(s) 
and  v(s)  corresponding  to  a  given  x  and  y. 

Either  complex  or  real  arithmetic  may  be  utilized  in  this 
program,  although  complex  arithmetic  offers  a  decided  advan¬ 
tage  in  convenience  and  simplicity.  However,  the  user 
should  be  aware  of  the  potential  dangers  in  unintentional 
exchange  of  sheets  of  a  Riemann  surface  when  investigating 
a  multi-valued  function.  More  will  be  said  about  this  later. 


Also,  although  FORTRAN  manuals  state  or  imply  that  the 
FUNCTIONS  DREAL  and  DIMAG  are  implemented  in  all  compilers, 
this  may  not  actually  be  the  case.  Accordingly  it  is  prudent 
to  add  the  following  statements  after  SUBROUTINE  VALUE. 

FUNCTION  DREAL (CPLX) 

REAL* 8  CPLX(2) 

DREAL  =  CPLX(l) 

RETURN 

END 

FUNCTION  DIMAG (CPLX) 

REAL *8  CPLX(2) 

DIMAG  =  CPLX ( 2 ) 

RETURN 

END 

SUBROUTINE  CURVE  is  also  a  user  supplied  routine  which 
calculates  the  x  and  y  values  for  each  consecutive  incrementa¬ 
tion  of  the  parameter  p  of  equation  (17),  thus  locating  a 
point  on  the  simple  parameterized  contour.  Equation  (18) 
provides  an  illustrative  and  highly  effective  example  of 
parabolic  form  which  has  been  employed  almost  exclusively 
during  this  investigative  work. 

SUBROUTINE  INCREP  formulates  the  increasing  sequence  of 
the  parameter  p  as  given  by  the  following  equation 

p  =  p  +  A  (20) 

Limited  analysis  with  other  than  equally  spaced  incrementa¬ 
tion  of  points  along  the  contour  of  integration  did  not  ex¬ 
hibit  any  enhancing  capabilities  and  was  not  continued. 

The  main  program  accomplishes  the  integration  of  equation 
(.14)  using  a  trapezoidal  summation  as  described  by  equation 
(16).  In  addition,  the  absolute  value  of  the  integrand  is 


summed  in  accordance  with  equation  (19).  Also,  the  number 
of  times  a  sign  changes  occurs  between  successive  evaluations 
of  the  addend  of  equation  (16)  is  recorded.  These  features 
afford  the  user  an  opportunity  to  monitor  potential  oscilla¬ 
tion  effects  upon  the  accuracy  of  the  results.  Clock  time 
and  the  number  of  calls  to  SUBROUTINE  VALUE  required  for 
performing  the  numerical  contour  integration  are  maintained 
as  measures  of  the  computational  efficiency  of  the  numerical 
inversion  procedure. 

Termination  of  numerical  integration  occurs  when  each 
of  N  successive  evaluations  of  the  addend  of  equation  (19) 
has  magnitude  less  than  a  specified  epsilon.  We  have  usually 
used  N  =  5.  This  is  a  very  stringent  requirement  and  could 
probably  be  relaxed  resulting  in  some  saving  of  computational 
time  and  with  but  minor  loss  of  accuracy. 

In  addition  to  the  termination  criterion  described  in 
the  previous  paragraph,  other  program  inputs  include  the 
increment  A  of  equation  (18),  appropriate  numerical  values 
of  any  constants  of  the  function  f(s),  the  starting  position 
A  of  equation  (18)  and  the  value  of  t. 

The  program  output  prints  the  starting  position  of  the 
integration  contour,  the  values  given  by  equations  (16)  and 
(19),  the  final  values  of  x  and  y  at  the  termination  of  inte¬ 
gration,  the  number  of  changes  of  sign  between  successive 
addends  in  equation  (16),  the  total  number  of  calls  made  to 
SUBROUTINE  VALUE  during  the  process,  and  the  clock  time  re¬ 
quired  for  numerical  integration. 
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As  an  illustration  consider  the  following  example  where 


f(s) 


s 


s2  * 


(21) 


A  suitable  form  for  SUBROUTINE  VALUE  is  as  follows: 

SUBROUTINE  VALUE  (X,  Y,  T,  AL,  U,  V,  MMMMM) 

C  THIS  IS  TABLE  ENTRY  NUMBER  09  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REAL* 8  ( A,B,D-H,0-Z) 

IMPLICIT  COMPLEX* 16  (C) 

MMMMM  -  MMMMM  +  1 
ZERO-O. D+0 
CS-DCMPLX(X.Y) 

C AL-DCMPLX ( AL , ZERO ) 

CT*DCMPLX( T , ZERO ) 

CDEN-CS  *  *2 +CAL  *  *  2 
C-CS/CDEN 
CEXP»CDEXP(CS*CT) 

C»C*CEXP 
U=DREAL(C) 

V-DIMAG(C) 

RETURN 
END 

If  the  contour  indicated  by  equation  (18)  were  chosen  and 
numerical  inversion  were  performed  upon  the  f(s)  given  by 
equation  (21)  for  the  case  where  A  *  0.3,  A  *  0.1,  ot  =  0.125, 

N  =  5  and  e  =  l.D-11,  the  evaluation  of  F(2tt)  would  yield  the 


(Note:  The  parameter  a  is 
denoted  by  the  FORTRAN 
variable  AL. ) 


following  output 


AA 

F(2ir ) 

|  G(2ir )  |  Xp  Yf 

LOSC 

0.3 

7.07107D-1 

1 . 32635D+0  -5.46  2.4 

4 

MMMMM 
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where  AA,  LOSC,  and  MMMMM  are  the  FORTRAN  program  names  for 
the  contour  starting  position,  the  number  of  sign  changes 
between  successive  addends,  and  the  total  number  of  calls 
made  to  SUBROUTINE  VALUE,  respectively.  Similar  results 
for  the  almost  100  functions  f(s)  tested  during  the  course 
of  this  research  are  listed  in  Appendix  B. 
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IV.  AUXILIARY  SUBROUTINES 


Chapter  III  briefly  mentioned  the  potential  dangers  in¬ 
herent  in  the  utilization  of  complex  arithmetic  within  SUB¬ 
ROUTINE  VALUE  when  the  function  f(s)  is  multi-valued. 

If  f(s)  has  branch  points  or  winding  points  it  is  impor¬ 
tant  to  assure  that  one  remains  on  the  same  branch  of  the 
function  as  one  proceeds  along  the  path  of  integration. 

Available  FORTRAN  operations  do  not  assure  this  and  it  is 
necessary  to  employ  specially  programmed  algorithms.  The 
present  discussion  is  limited  to  the  two  simple  cases,  rais¬ 
ing  of  a  complex  number  to  a  non-integer  power  and  taking 
the  logarithm  of  a  complex  number. 

The  essence  of  the  matter  lies  in  the  FORTRAN  operation 
DATAN2(y,x),  or  its  equivalents,  which  must  be  relied  upon 
to  provide  the  argument  of  a  complex  number  in  polar  form. 

This  operation  always  provides  a  result  in  the  range 
— tt  <  9  <_  tt,  so  that  if  Arg  (z)  actually  passes  through  the 
values  (2n+l)n,  for  integer  n,  the  FORTRAN  produced  value 
of  Arg(z)  experiences  a  jump  of  ±  2ir.  If  one  is  determining 
the  logarithm  of  z  or  a  non-integer  power  of  z,  shifting  of 
branch  will  take  place  unless  the  continuity  of  Arg(z)  is 
restored. 

SUBROUTINES  CPOWER  and  CLOG  provide  for  this  continuity 
by  monitoring  changes  in  Arg(z)  and  adding  or  subtracting 
2nif,  whenever  it  is  appropriate  to  do  so,  to  the  value  of 
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DATAN2(y,x) .  As  they  are  presenting  written,  each  provides 
for  up  to  five  separate  calls,  within  SUBROUTINE  VALUE,  to 
either  SUBROUTINE  CPOWER  or  SUBROUTINE  CLOG. 

These  subroutines  are  listed  below, 


SUBROUTINE  CPOWER  ( 11  N,  ZCUT,  P,  J ,  N C ALU 
IMPLICIT  REAL*8  (A-H,0-Z) 

COMPLEX  *  16  ZIN,  ZOUT 
DIMENSION  T0LD(5)  »  NCALL (5)  »  M D  I  (  5  ) 
IF(NCALLIJ) .EQ.l)  NPI(J>=0 
OH  =  0.0+0 
X  =  DREAL(ZIN) 

Y  =  DIM  AG( Z  IN  ) 

R  =  DSQRT  (  (X*X)  +  (Y*Y)  ) 

TR0N3  *  DATAN2(Y,X) 

IF( NCALL( J > .EO . 1 >  GO  TO  5 
PROD  =  ”rR0NG*T0LD <  J  ) 

IF(  PP0D.GT.0H.0R.  X.GT.OH)  GO  tq 
IF  (To  ONG.LT  .OH )  NPI(J)  =  NPI(J) 

--  -  N  pi <  Ji  =  NP  I  ( J  ) 


'  •  OHJ 


5 

+  2 
-  2 


IF(TRONG.G‘ 

RP  s  R**0 

EN  =  NP  I  (J  ) 

TRITE  »  tronG  ♦  ( EN*3.14159265358979324D+0) 
TP  =  TRI  T5*P 
X  =  RP*DC0S  CT P i 
Y  *  R P*D SI  N  (TP) 

ZOUT  =  DC  MP  LX(  X,  Y ) 

T  OLD  ( J  )  =  TRONG 
NCALL ( J)  =  NCALL( J) 

RETURN 
ENO 


♦  1 


SUBROUTINE  CLOG  (ZIN,  ZOUT,  J,  "CALL) 
IMPLICIT  RE  AL*8  (  A-H,  0-Z ) 

C OMR LE X* 16  ZIN,  ZOUT 

DIMENSION  T0L  0(5),  MCALLt  5 ) ,  NPI ( 5) 

I F( MCALL (J).EO.I)  NPIIJ)  *  0 
■OH  *  O.D+O 
X  *  OREAL(ZIN) 

Y  =  DIMAG(ZIN) 

R  *  DSQRT( ( X*X)  +  ( Y*Y ) ) 

TRONG  *  DAT  AN2 ( Y , X ) 

IF(MC4LL(J).EQ.  1)  GO  TO  5 
PROD  *  TRONG*TOLD(J ) 
IF(PROD.GT.OH.OR.X.GT.OH)  GO  tq  5 
I F ( TRONG. LT. OH)  NPI(J)  *  N«>I(J)  ♦  2 
IF(TRONG.GT.OH)  NPI(J)  *  NPIIJ)  -  2 
REAL  *  DLOG(R) 

EN  *  NPI  (J) 

TRITE  *  T  RONG  +  ( EN*3 . 14159265358979324D+0) 
ZOUT  *  DCMPLX(R£AL, TRITE) 

TOLD(J)  *  TRONG 

«CALL  (J  )  «  MCALL  (  J )  ♦  1 

RETURN 

END 
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If  either  algorithm  is  to  be  utilized  by  SUBROUTINE  VALUE, 
the  first  requirement  is  to  zero  the  applicable  integer 
array,  NCALL  or  MCALL,  within  the  main  program.  Then,  in 
SUBROUTINE  VALUE,  the  calls  to  SUBROUTINE  CPOWER  are  of 
the  form 

CALL  CPOWER(ZIN,ZOUT,P,  J, NCALL) 

where  ZIN  is  the  COMPLEX*16  operand,  ZOUT  is  the  COMPLEX+16 
output,  P  is  the  REAL*8  power,  and  J  is  the  INTEGER*4  index 
which  tells  whether  this  is  the  first,  second,...,  call  to 
SUBROUTINE  CPOWER  within  SUBROUTINE  VALUE.  The  relationship 
obtained  is 

ZOUT  =  (ZIN)P  (22) 

SUBROUTINE  CLOG  is  used  similarly.  Calls  made  from  SUBROUTINE 
VALUE  are  of  the  form 

CALL  CLOG( ZIN, ZOUT, J, MCALL) 

where 

ZOUT  =  Ln(ZIN)  (23) 

J  plays  the  same  role  as  in  SUBROUTINE  CPOWER.  There  is  no 
P.  It  should  be  noted  that  the  array  MCALL  replaces  NCALL. 
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THE  STEEPEST  DESCENT  CONTOUR 


The  theoretical  advantages  of  a  steepest  descent  contour 
have  been  mentioned  previously.  Consider  a  contour  in  the 
complex  plane  along  which  v=constant.  As  shown  in  Figure 
2,  let  the  £  axis  be  tangent  to  the  contour  with  the  n  axis 
perpendicular  to  the  £  axis  in  the  same  sense  that  y  is  perpen 
dicular  to  x. 


*  REALAOAIS 


Figure  2.  Steepest  Descent  Contour.  Contour  along 

which  v  is  constant.  Axis  ?  is  tangent  at 
point  P.  Axis  n  is  normal. 


Recall  equation  (11), 

g(s)  -  est  f(s)  =  u(s)  +  iv( s) 

By  the  Cauchy-Riemann  relationships 

3u  3v  3v  3u 
7T  "  Tn  '  "  TZ  *  3T? 

Along  a  v  »  constant  contour 


(ID 


(24) 


Such  a  v  =  c  =  constant  contour  divides  the  plane,  locally 

at  least,  into  two  regions:  with  v  >  c  and  R2  with  v  <  c. 

3  v 

Consequently  along  such  a  contour  is  one  signed,  and, 

dT) 

3  u 

therefore,  is  one  signed,  u  -►  0  along  the  v  =  c  contour, 
in  the  direction  of  positive  5.  Therefore,  there  is  no 
oscillation  in  u  and,  obviously  also,  there  is  none  in  v. 

Since  u  -*■  0  and  v  -*■  0  along  a  suitable  integration  contour, 
the  selection  of  c  =  0  assures  that  this  results. 

This  second  program  to  be  discussed  does  partially  follow 
a  path  of  steepest  descent.  First  it  integrates  vertically 
upward  from  a  chosen  point  on  the  real  axis  (along  a  Bromwich 
contour)  until  intersection  with  a  specified  v  =  0  contour 
occurs.  Then  it  employs  a  tracking  algorithm  to  obtain  suc¬ 
cessive  points  on  this  contour  along  which  it  then  integrates, 
moving  to  the  left. 

SUBROUTINE  VALUE  has  the  same  role  as  that  described  in 
Chapter  III,  namely,  it  calculates  the  u(s)  and  v(s)  in 
accordance  with  equation  (11)  corresponding  to  a  given  x  and 
y.  For  the  initial  part  of  the  algorithm,  the  integration 
along  the  Bromwich  contour,  the  value  of  x  is  constant,  and 
y  is  increased  in  increments  by  Ay.  Numerical  integration 
in  accordance  with  equation  (16)  is  conducted  just  as  it  was 
done  in  the  previous  case. 

On  the  real  axis,  obviously,  v  *  0.  With  an  initial 
incrementation  in  the  vertical  direction,  v  assumes  a  value 
other  than  zero.  The  sign  of  v  is  recorded  and  the  incrementa¬ 
tion  process  is  continued  until  the  sign  of  v  changes.  This 
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constitutes  a  contour  crossing.  The  procedure  continues 
until  a  predetermined  number  of  such  crossings  are  obtained, 
and  then  a  simple  linear  interpolation  is  performed  to  locate 
the  specific  ordinate  value  y  for  which  v  =  0.  Usually  one 
might  wish  to  use  the  first  v  =  0  contour  above  the  real 
axis,  but  the  second,  third,  or  other  such  contour  may  be 
selected. 

The  last  (y^»vi)  values  prior  to  the  desired  contour 

crossing  and  the  first  (y.,v.)  values  after  crossing  are 

J  J 

utilized  in  the  linear  interpolation  process  to  predict  YNEW, 
for  which  v  =  0.  SUBROUTINE  VALUE  provides  the  actual  VNEW 
value  corresponding  to  YNEW.  Depending  upon  the  sign  of 
VNEW,  the  appropriate  coordinate  pair  (3^*^)  or  (Yj.Vj)  is 
replaced  by  (YNEW, VNEW)  and  the  interpolation  is  repeated 
with  this  refinement.  Termination  of  the  interpolation  occurs 
when  YNEW  results  in  a  magnitude  of  VNEW  which  is  less  in 
magnitude  than  a  prescribed  tolerance  limit. 

Once  the  values  of  x,y,u,  and  v  have  been  established 
for  the  first  point  on  the  desired  v  =  0  contour  in  the  pre¬ 
vious  process,  the  second,  or  tracking  portion  of  the 
algorithm  commences. 

SUBROUTINE  XMARCH  is  called  to  provide  a  new  coordinate 
pair  (x,y)  in  the  negative  sense  of  the  real  axis  from  the 
known,  or  old,  coordinate  pair  (x0*y0)»  This  is  accomplished 
by  the  following  relationships 

x  =>  xQ  -  R  cos(0)  (26) 

y  ■  yQ  +  R  sin(0)  (27) 
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where  (R,9)  are  the  polar  coordinates  of  a  search  pattern 
measured  from  the  known  coordinates  (x0»y0)  on  the  v  ■  0 
contour.  The  angle  9  is  set  at  zero  radians  for  the  first 
call  to  SUBROUTINE  XMARCH.  SUBROUTINE  VALUE  is  then  called 
to  provide  the  value  of  v  at  this  point.  The  sign  of  v 
establishes  upon  which  side  of  the  v  =  0  contour  the  new  point 
lies. 

SUBROUTINE  ROTATE  is  then  called  to  increment  the  angle 
theta  in  an  appropriate  sense.  A  sweeping  search  is  then 
conducted  by  SUBROUTINES  ROTATE,  XMARCH,  AND  VALUE  until 
a  sign  change  is  detected  for  v  values  between  increments, 
thus,  indicating  a  crossing  of  the  v  =  0  contour. 

The  last  (0i»vi)  coordinate  pair  prior  to  a  crossing  and 

the  first  (B.jVj  coordinate  pair  after  a  crossing  are 
J  J 

utilized  in  a  linear  interpolation  scheme,  analogous  to  that 
of  the  first  program  segment,  to  refine  the  location  of  the 
new  point  on  the  v  »  0  contour. 

Once  the  second  point  on  the  v  *  0  contour  is  established, 
the  tracking  algorithm  may  be  enhanced  by  predicting  the 
third  point  on  the  contour  to  lie  at  the  same  angle  9  from 
the  second  point.  As  before,  the  sign  of  v  at  this  new  loca¬ 
tion  permits  a  rotational  search  to  be  conducted  in  the 
appropriate  sense.  For  small  values  of  the  radius  vector 

R  (R-.025  is  built  into  the  program,  but  may  be  changed  by  the 
user),  this  approximation  has  been  found,  in  general,  to  be  a 
reasonable  one. 

The  addends  of  equations  (16)  and  19)  are  evaluated  at 
each  step  in  the  marching  process  within  the  main  program 


and,  once  again,  termination  of  the  numerical  integration 
occurs  when  each  of  N  successive  addend  contributions  is 
less,  than  a  prescribed  epsilon. 

Required  inputs  to  the  main  program  are  the  index  speci¬ 
fying  the  desired  v  =  0  contour,  the  starting  position  on 
the  real  axis,  the  increment  of  the  vertical  march  along  the 
Bromwich  contour,  the  termination  criterion,  any  constants 
associated  with  the  function  f(s),  and  the  value  of  t.  These 
inputs  are  represented  by  the  program  variables  NREG,  AA, 
DELY,  EPS,  (AL,BE,...),  and  T,  respectively. 

The  simple  parameterized  contour  offers  one  decided  ad¬ 
vantage  over  the  steepest  descent  contour:  each  call  made  to 
SUBROUTINE  VALUE  is  productive,  in  the  former  case,  in  the 
sense  that  it  is  utilized  to  compute  a  contribution  of  the 
addend  of  equation  (16).  The  latter  algorithm  requires  calls 
to  SUBROUTINE  VALUE  to  perform  the  linear  interpolation  pro¬ 
cedures,  which  are  non-productive  in  terms  of  addend  computa¬ 
tion.  This  represents  increased  computer  time.  The  present 
investigation  has  shown  that  in  each  of  the  numerous  cases 
which  have  been  examined,  the  simple  parabolic  path  indicated 
by  equation  (18)  is  more  efficient  and  fully  as  accurate  as 
a  steepest  descent  contour. 
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VI .  SHANKS'  ACCELERATOR 


The  effectiveness  of  a  family  of  non-linear  sequence  to 
sequence  transformations  in  accelerating  the  convergence  of 
(some)  slowly  convergent  sequences  and  in  inducing  conver¬ 
gence  in  (some)  divergent  sequences  has  been  reported  by 
Shanks  [19].  If  {An>  (n=0,l,2, . . . )  is  a  sequence  of  numbers, 
we  form  a  sequence  of  sequences  (A^  n>  which,  for  convenience 
we  write  as  {A(k,n)>.  The  integer  k  indicates  the  "order" 
of  the  sequence,  with  Afl  =  A(0,n);  the  integer  n  indicates 
the  position  of  the  term  in  the  sequence.  The  rule  for  form¬ 
ing  the  sequence  {A(k+l,n)>  from  the  sequence  (A(k,n)}  is 

A(i,j)  =0  if  j  <  2i 


A(k+1  ,n  )  =  A<k?n.>.4(iL,n-2)  -.Alk.n-l)  (28) 
A(k,n)+A(k,n-2)-2A(k,n-l) 

n*2k+3, 2k+4, .... 

To  illustrate  the  power  of  this  computational  device  con¬ 
sider  the  application  of  equation  (28)  and  its  iterates  to 
the  first  ten  terms  of  the  alternating  series 


Ln(  2) 


1 

2 


1 

3 


1 

4 
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The  results  are 


A(l,n) 


1  WWIJ! 


A(2,n)  A(3,n) 


A(4,n) 


l. oooccocooo 
o.sooooooooo 

0.8233333333 
0. 5823333333 
0.7833333323 
0.  £i666666t>7 
.  C. 7595238095 

8  0.6345233095 

9  0.7456  349  20  6 
10  0.6456349206 


o.o 

00  0.0  0.0 
0I7000000000  0.0  0*0 

iMMi  CifHI?  SnuMai  0:0 

mm  tewtt 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


where  the  tenth  partial  sum,  A(0,10),  is  correct  to  only  one 
significant  figure.  However,  iterate  A(4,10)  is  correct  to 


eight  figures  since  Ln(2)  =  0.6931471806. 

In  this  thesis  research,  Shanks'  accelerator  has  been 
tested  extensively  in  conjunction  with  the  simple  parameter¬ 
ized  contours  of  integration  discussed  in  Chapters  II  and 
III. 

Modification  of  the  basic  algorithm  of  Chapter  III  con¬ 
sisted  of  the  elimination  from  the  main  program  of  the  ter¬ 
mination  criterion  for  the  numerical  integration  given  by 
equation  (16),  and,  instead,  performing  the  integration  for 
a  finite  number  of  terms.  These  addends  of  equation  (16) 
were  stored  in  a  linear  array  and,  subsequently,  transformed, 
within  the  main  program  in  accordance  with  equation  (28). 

The  accelerator  did  not  enhance  the  rate  of  convergence 
of  the  numerical  inversion  process  for  any  of  the  cases 
investigated.  In  fact,  if  the  number  of  terms  of  Aq,  which 
are  transformed,  becomes  sufficiently  large  so  that  the 
summation  of  the  addend  contributions  of  equation  (16) 
approaches  the  correct  result,  the  transformation  will 
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decrease  the  accuracy  of  the  result  as  the  denominator  of 
equation  (28)  becomes  small. 

| 

Appendix  A  contains  the  program  listing  for  the  simple 

I  I 

parameterized  contour  of  integration  in  conjunction  with 
Shanks'  accelerator. 

Shanks  [19]  discusses  transformations  of  higher  order. 
These  transformations  may  be  worthy  of  investigation. 


VII.  LOCATION  OF  FUNCTION  SINGULARITIES 


The  criteria  of  Chapter  II  for  the  distortion  of  the 
Bromwich  contour  include  the  requirement  that  all  of  the 
singularities  of  a  function  f(s)  be  to  the  left  of  the  con¬ 
tour.  Heretofore,  the  discussion  has  assumed  that  the  loca¬ 
tion  of  these  singularities  was  known.  Indeed,  this  is  not, 
in  fact,  always  the  case.  In  this  chapter  we  investigate  a 
method  of  determining  the  location  of  the  singularities  of 
f(s)  when  such  information  is  either  not  known  or  is  diffi¬ 
cult  to  obtain  analytically. 

To  see  how  one  can  proceed  without  knowing  where  the 
singularities  are,  consider  the  simple  case 


f  (s) 


~2  ”  2 
s  +  ot 


(30) 


for  the  particular  values  of  a  *  0.125  and  t  =  2ir . 

We  use  a  simple  parabolic  contour  given  by  equation  (13) 
with  B*1  and  with  A,  which  we  denote  by  the  symbol  xq  in  the 
discussion  to  emphasize  its  actual  significance,  being  given 
various  values.  If  we  plot  the  result  of  the  inversion 
versus  the  starting  position  x  ,  the  result  is  as  shown  in 
Figure  3. 

For  values  of  xq  less  than  zero,  the  result  is  0.00. 

For  xQ  in  the  range  from  about  0.2  to  2.0,  the  result  is 
0.707107.  For  xQ  larger  than  about  2.0,  the  result  is  near 
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For  x  in 
o 


0.7  but  is  sensitive  to  the  precise  value  of  xq. 
the  range  from  0.0  to  0.2,  the  result  exhibits  a  "spike" 
and  a  very  significant  change  from  0.0  to  0.707107. 


Figure  3. 


F(t)  vs  contour  starting  position  for 
cos(at).  Value  of  inverse  of  s 


for  t  a  2ir  as  a  function 
of  starting  abscissa  xQ. 


>lfe  can  explain  this  behavior  as  follows.  For  xQ  less 
than  zero,  the  singularities  are  on  the  wrong  side  of  the 
contour  and  we  simply  get  the  wrong  answer.  For  xQ  in  the 
neighborhood  of  0.0  to  0.2,  the  contour  passes  so  close  to 
the  singularity  that  the  trapezoidal  integration  scheme  is 
not  accurate  and  experiences  difficulty  in  getting  an 
accurate  value  either  for  the  right  answer,  0.707107  or  for 
the  wrong  answer  0.0.  For  xQ  in  the  range  from  0.2  to  2.0 
the  singularities  are  on  the  correct  side  of  the  contour  and 
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we  get  the  correct  answer.  For  xQ  larger  than  2.0,  there  are 
numerical  difficulties,  probably  due  to  oscillation  or  to 
the  effects  of  the  exp(st)  factor  for  s  having  large  real 
part. 

We  have  done  the  same  thing  with  many  other  transforms 
and  have  found  the  same  behavior.  Our  experience  leads  to 
this  suggestion  for  cases  where  the  user  does  not  know  the 
location  of  the  singularities.  Perform  the  calculation  for 
a  number  of  different  values  of  xQ  and  plot  the  results  as 
in  Figure  3.  The  correct  result  is  the  ordinate  of  the 
rightmost  plateau,  before  oscillations  or  other  numerical 
difficulties  have  had  a  chance  to  introduce  inaccuracies. 

In  the  case  illustrated  in  Figure  3,  there  were  only  two 
plateaus,  one  corresponding  to  all  the  singularities  lying 
on  the  wrong  side  of  the  contour  and  one  corresponding  to 
all  the  singularities  lying  on  the  correct  side  of  the  con¬ 
tour.  However,  with  more  singularities,  one  might  expect 
to  have  several  plateaus,  only  the  rightmost  of  which  repre¬ 
sents  the  correct  result. 
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VIII.  THE  HYPERBOLIC  TANGENT 


In  Chapter  II  it  was  stated  that  great  advantage  accrues 
from  distorting  the  Bromwich  contour  as  long  as  the  process 
of  distortion  does  not  cause  the  contour  to  cross  over  a 
singularity  of  f(s).  In  this  chapter  we  investigate  a  typi¬ 
cal  case  in  which  it  is  impossible  to  satisfy  this  condition. 

The  function 


f (s)  =  i  tanh(~)  (31) 

has  an  infinite  sequence  of  simple  poles  on  the  imaginary 
axis  at  points  s  =  ±(2n+l)^  .  By  other  methods  it  may  be 
shown  that  this  function  is  the  transform  of  the  square  wave 
function  shown  in  Figure  4.  In  this  chapter  we  investigate 
the  extent  to  which  our  present  methodology  is  capable  of 
obtaining  this  square  wave. 
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Figure  4.  Square  wave  function. 

The  convergence  of  the  numerical  integration  has  been 
assured  by  bending  the  contour  to  the  left  in  all  other  cases 
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reported  in  this  thesis.  Accordingly,  we  do  so  in  the  pres¬ 
ent  instance.  This  implies  that  we  must  pass  between  two 
poles  as  the  contour  bends  to  the  left  and  all  of  the  in¬ 
finite  number  of  poles  above  the  point  where  the  contour 
crosses  the  imaginary  axis  thus  fall  on  the  wrong  side  of 
the  contour  and  are  not  included  on  the  left.  Thus,  in 
effect,  our  procedure  is  finding  the  inverse  Laplace  trans¬ 
form  of  a  substitute  function  which  has  only  a  finite  number 

of  poles  coinciding  in  location  with  those  of  the  given 

\ 

function  f(s)  and  having  the  same  residues.  There  is  some 
reason  to  hope  that  the  inverse  of  the  substitute  function 
will  be  essentially  similar  to  that  of  the  given  function, 
at  least  for  sufficiently  small  values  of  t. 

Indeed  we  do  find  this  to  be  the  case.  SUBROUTINE  CURVE 
of  the  algorithm  of  Chapter  II  was  modified  to  formulate 
a  simple  parameterized  curve,  which  utilized  the  Bromwich 
contour  for  numerical  integration  in  its  initial  segment, 
starting  from  a  position  s-y  on  the  positive  real  axis  and 
continuing  until  the  elevation  of  a  finite  number  of  poles 
is  reached.  The  second  segment  of  the  simple  parameterized 
curve  which  was  used  consists  of  the  parabola  of  equation 
(18).  -This  parabolic  arc  allowed  the  integration  contour 
to  cross  the  imaginary  axis  between  adjacent  poles  of  the 
function  f(s).  Such  a  contour  is  indicated  in  Figure  5. 

The  termination  criterion  which  was  used  in  this  case 
is  identical  to  that  described  in  Chapter  III.  All  other 
program  input  values  remained  unaltered  from  the  form  of  their 
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Figure  5.  Contour  of  integration  for  the  hyperbolic 
tangent..  Special  contour  for  inversion  of 
f(s)  =  T)  tanh(as/2).  It  starts  at  xq, 
rises  vertically,  and  then  follows  a  para¬ 
bolic  path  between  two  poles. 

previous  presentation  with  tl.c  exception  that  the  integer 
JJJJJ,  which  selected  the  number  of  function  poles  to  be 
encompassed  by  the  contour  was  added.  The  program  listing 
appears  in  Appendix  A. 

The  approximate  solutions  F(t)  obtained  by  numerical 
inversion  of  the  substitute  function  for  the  f(s)  given 
by  equation  (31)  are  shown  graphically  below  in  Figures  6-9. 

Figure  6  is  a  reproduction  of  five  periods  of  the  square 
wave  of  Figure  4  with  a  =  10.0  and  the  first  100  poles  of 
f(s)  above  the  real  axis  to  the  left  of  the  contour  of  in¬ 
tegration.  Figure  7  is  a  similar  picture  of  the  same  F(t) 
with  a  =  5.0  and  the  first  150  poles  of  f(s)  above  the  real 
axis  encompassed  by  the  contour,  with  only  one  period  of  the 
wave  form.  Figures  8  and  9  are  analogous  to  Figures  6  and  7, 
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Figure  6. 


Numerical  approximation  of  five  cycles  of  the 
square  wave  function  with  100  poles  to  the 
left  of  the  contour. 
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Figure  7.  Numerical  approximation  of  one  cycle  of  the 
square  wave  function  with  150  poles  to  the 
left  of  the  contour. 


respectively ;  however,  only  the  first  seven  poles  of  f(s) 
above  the  real  axis  were  included  on  the  correct  side  of  the 
contour  of  integration. 

The  reader  should  be  aware  that  the  increment  of  t  used 
in  the  above  series  of  graphs  differed  in  each  case.  The 
inaccuracies  shown  in  Figures  6-9  are  partially  due  to  the 
assumptions  inherent  within  the  employment  of  our  algorithm, 
and,  also,  partially  due  to  the  fact  that  the  plotter  simply 
connects  points  sequentially.  This  can  also  be  seen  in  the 
slope  of  the  graphs  at  t  =  a,  2a,  etc.,  which  is  a  function 
of  the  particular  At  employed  and  not  an  inaccuracy  of  the 
algorithm. 

Efforts  to  accomplish  numerical  inversion  of  the  f(s) 
of  equation  (31)  using  the  steepest  descent  contour  were 
unsuccessful.  The  difficulty  in  using  this  contour  appeared 
to  be  related  to  an  adverse  influence  from  the  proximity  of 
the  function  poles  to  the  v  *  0  contours  near  the  imaginary 
axis.  The  algorithm  was  overtaxed  in  this  vicinity  and 
could  not  track  the  contour;  consequently,  further  investi¬ 
gation  was  not  attempted. 

In  summary,  we  conclude  that  the  effect  of  substituting 
a  function  which  has  only  a  finite  number  of  poles  which  are 
coincident  in  location  with  the  first  n  poles  of  a  function 
f(s),  which  possesses  an  infinite  number  of  poles  spaced  incre¬ 
mentally  along  the  imaginary  axis,  produces  an  inverse 
approximation  to  the  inverse  of  the  given  function  f(s). 

The  accuracy  obtained  with  such  a  substitution  is  increased 
as  more  of  the  function  singularities  are  encompassed  by  the 
distorted  contour. 

. . . . . m . m . . . . . 


IX.  FACTORS  WHICH  INFLUENCE  ACCURACY 


Heretofore,  our  discussion  has  included  the  treatment  of 
contour  integration,  utilizing  a  simple  parameterized  curve 
for  that  purpose,  from  two  perspectives.  The  first  of  these 
methods,  discussed  in  Chapter  III,  incorporated  a  termina¬ 
tion  criterion  which  required  a  prescribed  number,  N,  of 
successive  integrand  contributions  each  to  have  smaller  mag¬ 
nitude  than  a  specified  epsilon.  The  second  of  these  methods, 
discussed  in  Chapter  VI,  performed  numerical  integration  for 
a  finite  number  of  points  spaced  along  the  contour  with  sub¬ 
sequent  employment  of  Shanks'  accelerator  in  an  endeavor  to 
enhance  convergence. 

The  first  of  these  methods,  we  have  found,  was  highly 
successful  when  employed  with  the  proper  point  spacing,  and 
a  very  tight  termination  criterion,  namely,  a  small  epsilon 
and  large  N.  The  second  method  has  not,  in  our  investiga¬ 
tions,  been  found  to  be  successful,  and  need  not  be  discussed 
further.  The  issue  then  becomes  one  of  selecting  the  proper 
combination  of  step-size,  epsilon,  and  N  to  provide  an 
accurate  and  economical  result. 

In  order  to  make  such  an  appropriate  selection  of  these 
parameters,  it  is  necessary  to  examine  the  factors  which  in¬ 
fluence  accuracy.  These  factors  include:  the  step-size  along 
the  contour,  the  termination  criteria,  and  the  oscillation. 


Consider  a  simple  case,  which  we  have  chosen,  along  a 
path  for  which  the  oscillation  is  insignificant.  The  accu¬ 
racy  in  such  a  case  is  influenced,  if  the  termination 
criteria  are  sufficiently  stringent,  only  by  the  step-size 
along  the  contour  of  integration.  Thus  dealing  with  the 
transform 

f(s)  =  y—  "2  (32) 

,  s'  +  a 

whose  inverse  is  the  cosine  function,  yields  a  plot  of  the 
logarithm  of  the  percent  error  of  F(t)  versus  the  step-size, 

A,  as  shown  in  Figure  10.  In  this  case,  the  contour  of  in¬ 
tegration  was  the  simple  parameterized  curve  of  equation 
(18)  with  A  =  0.3,  B  =  1.0,  a  =  0.125  and  t  =  2tt. 

The  sequences  of  points  in  this  plot  are  all  obtained  by 
maintaining  N  constant  (N=3)  and  plotting  log10  Percent  error 
versus  A  with  epsilon  as  a  parameter. 

Four  such  parameter ic  plots  appear  in  the  figure  below 
corresponding  to  the  values  of  epsilon  (EPS)  as  shown  in 
the  legend. 

Clearly  the  accuracy  of  the  results  may  be  seen  to  behave, 
in  a  general  sense,  in  the  manner  shown  in  Figure  11,  where 

EPSj^  <  eps2  <  eps3  <  eps4. 

When  the  termination  criteria  are  tight,  namely,  when  EPS 
is  small  and  N  is,  at  least,  greater  than  one,  the  percent 
error  is  a  narrowly  banded,  roughly  linear  function  of  the 
integration  step  size,  A,  over  a  range  from  about  A  =  0.02 
to  A  »  0.2.  This  corresponds  to  EPSj^  in  Figure  11. 
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Figure  10.  Logarithm  of  percent  error  versus  step  size 
for  the  cosine  function. 


If  the  termination  criterion  EPS  is  looser,  as  evidenced 
in  other  cases,  shown  schematically  in  Figure  11,  the  accu¬ 
racy  of  the  result  cannot  be  improved  indefinitely  by  decreas¬ 
ing,  A.  Rather,  there  is  a  point,  which  differs  in  each  case, 
beyond  which  the  accuracy  is  decreased  as  the  step  size  is 
decreased.  The  behavior  is  complicated  by  the  fact  that  the 
effects  of  changes  in  step-size  and  changes  in  termination 
criteria  are  not  themselves  disjoint.  Making  the  step-size 
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Schematic  representation  of  accuracy  factors 
Behavior  of  the  logarithm  of  percent  error 
versus  step  size  A  as  a  function  of  the  ter¬ 
mination  criterion  EPS. 


smaller  results  in  meeting  the  termination  criteria  at  an 
earlier  point  along  the  path  of  integration  and  thus  may 
actually  reduce  overall  accuracy. 

One  way  of  dealing  with  this  matter,  of  course,  is  to 


revise  the  termination  criteria  so  as  to  remove  this  inter 


dependence.  We  simply  have  not  experimented  with  using 
alternate  forms  for  the  termination  algorithm.  Instead,  we 
suggest  the  following  viewpoint  to  the  user  who  wishes  to 


assure  that  he  obtains  results  with  a  specified  accuracy. 

He  should  first  of  all  employ  a  very  strict  termination 


criterion  and  use  a  succession  of  rather  small  values  of 


step  size,  and  should  also  vary  the  starting  point  on  the 
contour.  In  this  way  he  should  be  able  to  establish  an 
evaluation  having  even  greater  accuracy  than  he  requires. 
Then,  selecting  what  seems  to  be  a  good  starting  point  and 
keeping  to  it,  he  should  increase  the  step  size  until  the 
result  is  no  more  accurate  than  he  requires.  He  should 
observe  how  much  saving  in  computer  time  (or  in  the  number 
of  calls  to  VALUE)  this  affords  and  should  not  increase  step 
size  beyond  a  reasonable  point.  Then  he  should  loosen  the 
termination  criterion  progressively,  measuring  the  saving 
in  computer  time  versus  the  possible  loss  of  accuracy, 
stopping  short  of  the  point  where  he  cannot  rely  on  obtain¬ 
ing  the  required  accuracy. 

In  many  cases  where  numerical  integration  is  employed 
to  obtain  a  result  whose  correct  value  is  Y,  the  numerically 
produced  approximation  y  approximately  satisfies  a  relation 

y  *  Y  +  aAm  (33) 

where  a  and  m  are  unknown  constants.  The  exponent  m  may 
frequently  be  established,  once  and  for  all,  for  the  type 
of  integration  being  used,  and,  for  individual  integrands, 
the  correct  value,  Y,  from  two  evaluations,  y^  obtained  with 
A  *  A^  and  y^  obtained  with  A  =  Ag,  by  using  the  extrapola¬ 
tion  formula 

Y  '  sps;  IV2  -  ViJ  <34> 
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We  have  investigated  this  method  of  obtaining  improved  accu¬ 
racy  for  our  inversions  and  have  found  that  with  the  inter¬ 
dependence  between  termination  criteria  and  step  size,  the 
numbers  a  and  m  were  both  rather  large  and  difficult  to 
establish  -  in  other  words,  the  extrapolation  was  not  success¬ 
ful  in  producing  improved  values.  Our  conclusion  was  that 
it  was  more  efficient  from  a  computational  point  of  view 
simply  to  use  a  sufficiently  small  step  size  and  an  appro¬ 
priately  matched  set  of  termination  criteria  so  as  to  be 
able  to  obtain  an  accurate  answer  without  interpolation. 

However,  a  suggestion  for  further  study  of  this  matter 
is  made  at  the  end  of  the  next  chapter. 
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RESULTS.  CONCLUSIONS.  AND  RECOMMENDATIONS 


This  thesis  has  demonstrated  that  numerical  integration 
along  a  suitable  contour  in  the  complex  plane,  thus  implement¬ 
ing  the  complex  inversion  integral  formula,  is  an  effective 
way  of  obtaining  numerical  values  for  the  inverse  Laplace 
transformation  of  a  given  function. 

The  success  of  the  procedure,  i.e.,  its  efficiency,  avoid¬ 
ance  of  inaccuracy  due  to  oscillation,  and  its  termination 
after  including  only  a  reasonable  number  of  addends  in  the 
numerical  sum,  is  related  to  the  fact  that  the  contours  chosen 
bend  to  the  left  as  they  rise,  thus  taking  advantage  of  the 
exp(st)  factor  in  the  integrand  of  the  inversion  formula. 

Although  there  are  some  theoretical  advantages  to  what 
we  have  called  the  "steepest  descent"  contour  which  follows 
a  v*0  path,  the  algorithm  which  permits  following  such  a  path 
requires  repeated  evaluations  of  the  integrand,  many  of  which 
are  not  actually  used  in  forming  the  summands  for  the  numeri¬ 
cal  summation.  We  have  found  by  investigating  almost  100 
cases  for  which  known  inverses  are  available  in  analytic  form, 
that  it  is  more  efficient  to  use  a  simpler  contour  than  the 
steepest  descent  path.  By  using  what  we  have  called  a  simple 
parameterized  curve,  we  have  devised  a  procedure  which  makes 
use  of  every  evaluation  of  the  integrand.  We  have  found  that 
the  penalties  of  not  following  a  steepest  descent  curve, 


namely  an  increase  in  oscillation  of  the  integrand  along  the 
path,  may  be  made  quite  tolerable  by  using  any  of  a  number 
of  simple  curves.  In  particular,  most  of  our  ^successful  work 
has  been  with  a  simple  parabola.  Our  algorithm,  however, 
includes  the  output  of  information  which  can  serve  to  alert 
the  user  to  the  possibility  of  inaccuracy  due  to  oscillation 
and  thus  suggest  to  him  that  he  might  select  a  different  form 
for  the  path. 

We  have  even  been  successful  in  obtaining  the  inverse  of 
a  function  which  has  an  infinity  of  poles  spaced  along  the 
imaginary  axis.  In  doing  so  we  violated  the  principle 
that  all  poles  of  the  transform  must  be  to  the  left  of  the 
contour.  Nevertheless  the  results  are  quite  satisfactory. 

Our  success  and  accuracy  have  been  so  gratifying  that 
we  venture  to  suggest  that  our  method  may  prove  to  be  an 
efficient  alternate  method  for  the  evalution  of  exotic  func¬ 
tions  for  which  other  methods  are  slowly  convergent  or  in¬ 
volve  series  the  coefficients  in  which  are  difficult  to 
obtain.  For  example,  case  92  in  Appendix  B  shows  the  suc¬ 
cessful  evaluation  of  the  rather  uncommon  Struve  function. 

The  original  reason  for  attempting  to  employ  numerical 
integration  in  the  complex  plane  as  a  means  of  inverting  a 
Laplace  transformation  was  that  alternate  methods  employed 
by  Hiep  and  Zargary  in  conjugated  heat  transfer  problems 
were  simply  not  accurate  enough.  They  led  to  physically 
impossible  results  with  some  temperatures  in  the  media  below 
sink  temperature  and  ochers  above  source  temperature.  . u 


this  thesis  we  have  not  restudied  these  heat  transfer  prob¬ 
lems  and  recommend  that  this  be  dore  using  our  methods  of 
inversion.  There  is  some  reason  to  believe  that  the  diffi¬ 
culty  encountered  by  Hiep  and  Zargary  will  not  be  encountered 
if  our  method  is  used,  but  we  are  not  prepared  to  substan¬ 
tiate  this  claim  at  this  time. 

Also  it  might  be  of  use  if  the  effects  of  varying  the 
available  parameters  (shape  of  curve,  spacing  of  points, 
starting  point  of  curves,  termination  criteria,  etc.)  were 
to  be  investigated  further  so  that  a  user  would  be  better 
prepared  to  deal  with  indications  of  inaccurate  inversion 
or  inordinate  requirements  for  computer  time.  Our  own 
numerical  experiments  were  invariably  so  happily  successful 
that  we  have  not  encountered  need  for  such  information. 

For  functions  f(s)  for  which  it  may  be  difficult  to 
locate  the  singularities,  we  have  shown  that  varying  the 
starting  point  of  a  simple  parameterized  contour  permits 
selecting  an  optimum  contour  in  the  sense  that  one  can  be 
assured  that  all  singularities  are  to  the  left  of  the 
curve  and  also  that  the  curve  does  not  reach  far  enough  to 
the  right  to  impose  numerical  difficulties  with  large  posi¬ 
tive  exponentials. 

At  the  end  of  the  preceding  chapter  we  indicated  that 
we  did  not  find  it  profitable  to  employ  extrapolation  as  a 
means  of  obtaining  improved  accuracy.  However,  this  is 
probably  worth  looking  at  again,  and  our  suggestions  for 
doing  so  are  as  follows.  First  employ  a  termination  criterion 


50 


which  is  disjoint  from  the  integration  step  size  in  the 
sense  that  the  location,  along  the  path  of  integration,  of 
the  point  where  termination  occurs,  is  independent  of  step 
size.  One  way  of  doing  this  is  to  employ  an  epsilon  in  the 
termination  criterion  which  is  itself  the  result  of  multiply¬ 
ing  a  fixed,  input  epsilon  by  the  step  size.  Then  the 
extrapolation  has  a  chance  of  being  successful.  So  as  to 
maintain  optimum  computational  efficiency,  one  should  use 
what  is  called  Richardson  extrapolation  in  which  the  evalua¬ 
tion  points  for  the  larger  step  size  are  themselves  included 
among  those  for  the  smaller  step  size. 
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APPENDIX  A 


LISTING  OF  COMPUTER  PROGRAMS  AND  SUBROUTINES 

The  computer  programs  and  principal  subroutines  that  we 
have  prepared  and  tested  during  our  investigation  are  all 
listed  within  this  appendix.  Section  A-l  contains  the  user's 
instructions  glossary,  flowchart,  and  program  listing  for  the 
simple  parameterized  curve.  'Section  A-2  is  a  similar  treat¬ 
ment  for  the  steepest  descent  contour.  Section  A-3  contains 
the  user's  instructions  and  program  listing  for  the  simple 
parameterized  curve  procedure  adapted  for  use  with  Shanks' 
accelerator.  Section  A-4  contains  the  same  information  for 
adaptation  of  the  simple  parameterized  curve  for  use  with 
the  hyperbolic  tangent  function  described  in  Chapter  VIII. 

It  is  hoped  that  sufficient  detail  has  been  provided  to 
enable  a  user  to  adapt  one  of  these  programs  Co  his  purpose 
in  an  efficient  manner.  However,  if  additional  insight  is 
required,  it  may  be  necessary  to  refer  to  the  chapter  of 
this  thesis  in  which  the  algorithm  is  developed.  This  is 
particularly  true  in  the  cases  of  Sections  A-3  and  A-4;  the 
development  in  these  sections  has  not  been  repeated  where  it 
is  equivalent  to  that  of  Section  A-l. 

Functions  DREAL  and  DIMAG  and  SUBROUTINES  CPOWER  and 
CLOG,  if  they  are  required,  are  listed  within  Chapter  III, 
as  is  a  typical  example  of  SUBROUTINE  VALUE.  Additionally, 
all  functions  f(s)  which  have  been  investigated  during  this 
research  have  their  applicable  SUBROUTINE  VALUE  listed  within 
Appendix  B. 


SECTION  A-l :  THE  SIMPLE  PARAMETERIZED  CURVE 

USER  INSTRUCTIONS 

This  is  the  FORTRAN  IV  program  for  the  numerical  inver¬ 
sion  of  the  Laplace  transformation  of  a  function  f(s),  which 
performs  contour  integration  along  a  distorted  Bromwich  con¬ 
tour  in  the  form  of  a  simple  parameterized  curve.  The  fol¬ 
lowing  instructions  are  provided  to  assist  the  user  in 
familiarizing  himself  with  the  program  so  that  he  can  adapt 
it  to  his  particular  requirements. 

1.  SUBROUTINE  CURVE  is  a  user  supplied  group  of  instructions 
which  formulates  the  numerical  integration  contour.  The 
program  listing  presently  contains  a  simple  parabola,  which 
may  be  utilized  directly  if  desired. 

2.  SUBROUTINE  VALUE  is  a  user  supplied  subroutine  which 
calculates  the  real  and  imaginary  components  for  g(s)  = 
exp(st)*f(s)  =  u(s)  +  iv(s).  For  examples  of  preparation 
of  this  subroutine  the  user  is  referred  to  the  test  cases 
of  Appendix  B. 

3.  The  mandatory  input  variables  to  the  main  program  are 
as  follows: 

(a)  AA  -  REAL*8  starting  position  (S=AA)  on  the  real 
axis  of  the  contour  of  integration.  (Also  denoted 
by  A  in  equation  (18)  and  xQ  in  Chapter  VII.) 

(b)  AL  -  REAL*8  value  of  any  constant  associated  with 
the  function  f(s).  (Note:  If  the  function  f(s) 
involves  more  than  one  constant,  the  calling  state¬ 
ments  to  SUBROUTINE  VALUE  must  be  modified  accord¬ 
ingly. 


(c)  DELP-REAL*8  increment  of  the  p-parameter  leading  to 
spacing  of  points  along  the  contour  of  integration. 

(d)  T  -  REAL*8  time  for  which  F(t)  is  desired. 

4.  Input  parameters  which  may  be  altered  by  the  user,  at 
his  discretion,  are  as  follows: 

(a)  EPS  -  REAL*8  tolerance  parameter  for  the  termina¬ 
tion  of  numerical  integration.  The  program  listing 
presently  contains  a  value  of  l.D-11. 

(b)  NUMBER  -  INTEGER*4  number  of  successive  addend  con¬ 
tributions  less  than  the  specified  value  of  EPS 
for  which  numerical  integration  is  terminated. 

5.  The  output  variables  which  are  printed  as  output  by  pro¬ 
gram  are  defined  below: 

(a)  AA  starting  position  on  the  real  axis 

(b)  SUM  the  value  of  F(t)  obtained  by  numerical 

inversion  of  f(s).  (cf.  equation  (16)) 

(c)  SUMA  the  absolute  value  of  the  summation  of  the 

the  addend  contributions  of  the  numerical 
integration,  (cf.  equation  (14)) 

(d)  X  the  final  value  of  x  on  the  contour  at  which 

integration  was  terminated 

(e)  Y  the  final  value  of  y  on  the  contour  at  which 

integration  was  terminated 

(f)  LOSC  the  number  of  sign  changes  between  successive 

addends  encountered  during  numerical  integra¬ 
tion 

(g)  MMMMM  the  total  number  of  calls  made  to  SUBROUTINE 

VALUE  during  the  numerical  integration  process. 
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6.  The  output  of  these  variables  is  as  follows: 

AA,  SUM,  SUMA,  X,  Y,  LOSC,  MMMMM 

7.  The  time  (seconds)  required  to  perform  the  numerical  inte¬ 
gration  is  printed  on  a  line  immediately  preceding  the  items 
listed  above. 

8.  As  previously  described,  the  program  will  perform  numeri¬ 
cal  integration  along  one  simple  parameterized  curve.  If  the 
user  desires,  the  program  may  be  modified  to  "sweep"  the  con¬ 
tour  of  integration  over  a  range  of  starting  values,  AA,  (cf. 
Chapter  VII).  This  may  be  quickly  accomplished  in  the  follow¬ 
ing  manner. 

(a)  Replace  the  value  of  the  integer  KSTART  from  1  to 
the  number  of  contour  integrations  desired. 

(b)  Replace  the  REAL*8  value  of  AA  with  the  following 
statement 

AA  =  XSTART+AK*XINCRM 

where 

XSTART  *  the  smallest  starting  value  of  the 
contour  of  integration  the  user 
desires 

XINCRM  *  the  Ax  between  starting  positions 

of  successive  contours  of  integration 
must  also  be  supplied^ 

Thus,  the  program  will  perform  numerical  inversion  along 
KSTART  contours  of  integration  which  differ  only  in  their 
starting  positions  along  the  x  axis.  The  user  may  then 
observe  increases  or  jumps  in  the  value  of  the  output  variable 
SUM  between  successive  integrations,  enabling  him  to  ascertain 
the  locations  of  function  poles.  The  utility  of  this  proce¬ 
dure  is  discussed  in  depth  in  Chapter  VII. 
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9.  If  SUBROUTINES  CPOWER  and/or  CLOG  are  called  by  SUBROU¬ 
TINE  VALUE  in  conjunction  with  use  by  this  program,  NCALL 
and/or  MCALL  must  be  dimensioned  and  initialized  to  zero 
wlchin  the  main  program  and  added  to  the  calling  statements 
to  SUBROUTINE  VALUE.  Additionally,  any  powers,  other  than 
integer  values,  must  be  entered  within  the  main  program  and 
passed  within  the  calling  statements  to  SUBROUTINE  VALUE. 


“ - 5 . . . . . *=■= 

*1 

GLOSSARY  OF  VARIABLE  NAMES 

VARIABLE 

DESCRIPTION 

LEGEND 

AA 

starting  position  on  the  real  axis 
of  the  simple  parameterized  curve 

2 

ADD 

•  i 

computed  value  of  the  addend  of  the 
numerical  integration 

1 

ADDA 

l  »• 

: 

i 

absolute  value  of  the  addend  of  the 
numerical  integration 

1 

ADDDLD 

saved  value  of  the  previous  addend 
contribution 

1 

l  AL 

constant  associated  with  f(s) 

2 

DELP 

increment  of  point  spacing  along  the 
simple  parameterized  curve 

2 

EPS 

tolerance  parameter  for  termination 
of  the  numerical  integration  process 

4 

KOUNT 

integer  counter  utilized  to  terminate 
the  numerical  integration  process 

1 

,  KSTART 

I 

integer  counter  for  the  number  of  times 
the  numerical  integration  is  desired 
to  be  performed,  using  different  start¬ 
ing  positions 

4 

LOSC 

' 

integer  counter  utilized  to  record  the 
number  of  sign  changes  between  succes¬ 
sive  addends 

3 

MMMMM 

integer  counter  to  record  the  total 
number  of  calls  made  to  SUBROUTINE  VALUE 
during  the  numerical  integration 

3 

NN 

integer  counter  to  record  the  total 
number  of  addend  contributions 

1 

NUMBER 

ti 

integer  input  of  the  desired  number  of 
successive  addend  contributions  less 
than  a  specified  epsilon  for  termination 

4 

|  p 

parameter  of  the  simple  parameterized 
curve 

PROD 

1 

1 

flag  utilized  to  detect  oscillation  of 
the  integrand 

1 
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SUM  total  addend  contribution  for  the  3 

numerical  integration  process 

SUMA  absolute  value  of  the  total  addend  3 

contribution  for  the  numerical 
integration  process 

T  time  for  which  F(t)  is  to  be  evaluated  2 

TEST  flag  for  termination  of  the  numerical  1 

integration  process 

U  real  part  of  the  function  g(s)  1 

UO  previous  value  of  U  1 

V  imaginary  part  of  the  function  g(s)  1 

VO  previous  value  of  V  1 

X  real  part  of  s  1 

XO  previous  value  of  X  1 

Y  imaginary  part  of  s  1 

YO  previous  value  of  Y  1 


LEGEND 

1.  No  action  required  by  the  user 

2.  Mandatory  user  input 

3.  Appears  as  program  output 

4.  May  be  altered  by  the  user  at  his  discretion 


ESTABLISH  THE 
SYSTEM  LIMITS 
AND  PARAMETERS 
INITIALIZE 
SET  THE  TIMER 


CALL 

CURVE  AND  VALUE 
OBTAIN 
INITIAL 
U  AND  V 


DO  400 

INCREMENT  THE 
STARTING  POINT 
COMMENCE 
NUMERICAL 
INTEGRATION 


* 


IF  NN 


ESTABLISH 
THE  SIGN  OF 
THE  INTEGRAND 


COMPUTE  THE 
PRODUCT  OF  THE 
OLD  AND  NEW 
INTEGRANDS 


PRODUCT 


LOSC _ 

THE  COUNTER 
LT.O.DfO  FOR  THE 

> - - - >  NUMBER  OF 

SIGN  CHANGES 
IS  INCREMENTED 


STOP  THE 
TIMER 


* 


WRITE 

THE 

OUTPUT 


t 


400 

CONTINUE 


THIS.  IS.  THE  .PROGRAM  JC.. PE  W ORM.NUFERI  CA^  INVf  RS ICN.  CF.I FE.LAPLACJ 


TRANSFORMATION  OF  A  FUNCTION  USING  THE 
IMPLICIT  REAL*8  (A-H.C-Z) 

CALL  ERRSET  ( 207, 256,- 1. I,  i ,  2 OS) 

NUMBER  »  5 - - - — 

aJd+omdarsini i.d+oi 

AL  -  1.250-1 - 


[MPLE  PARAMETERIZED  CURVE 

LEGENC 

- -  c. 

- -  4 

_ _  6 


• 

•  THE  OUTER  i.00P  MAY  BE  UTILIZED  TC  SWEEP  THE  CCNTCUR  STARTING 
•POSITION  OVER  A  DESIRED  RANGE  AS  DESCRIBED  IN  ThE  LSER  INST. 

• 

00  400  KKK-l.KSTART 

ADO  *  0.0*0 

ADDA  «  C.0+0 

SUM  «  0.0*0 

SUMA  »  C.0+0 

LOSC  *  0 

NPMMM  »  0 

AK  «  K.KK  _ 

AA  -  3.D-1 - - 

T  -  TWOPI - «.  ’Z 

P  ■  0.0+0 
NN  ■  0 
KCUNT  ■  0 

EPS  -  1.0- U _ _ 


1 .0-L 


a 


OELP  .  —  - 
CALL  CURVE  ( AA , P « XO*  Y C) 

CALL  VALUE  (X3,Y3 ,T,AL,UO, VO.MMMMM) 

TIME  «  0.0+0 
CALL  SETIME 

s  $ 

*  THE  FOLLOWING  LOOP  PERFORMS  THE  NUMERICAL  INTEGRATION  * 

NN  ■  NN+l 

CALL  INCPEP  ( P t  OELP) 

CALL  CURVE  (AA,P,X,  Y I 

CALL  VALUE  (X,Y ,f  ,AL,U,V  ,MMMMM ) 

ADO  -  ((V+V0J*1X-X0)+(U+U0>*IY-Y0J)/TW0PI 

IF  (NN.EO. II  A003LJ«A00 

PROD  -  ADD* AOQOLD 

IF  I PROD.LT.G.D+UI  lOSC»13SC+1 

ACDA  «  CABS (ADD  I 

SUM  »  SUM+AOD 

SUMA  «  SLMA+ADDA 

*********»******•***********»****»*»****•****************+*******• 

*  * 

*  THIS  IS  THE  TERMINATION  CRITERION  SEGMENT  OF  THE  ALGORITHM  * 

TEST  ■  DABS ( ADOI 

IF  (TEST .LT  .EPS  )  KOUNT*KCUNT+ 1 

IF  (TEST. G£. EPS!  KOUNT«0 

CALL  GET  IM  E  ( is  TJ 

EL  -  DFLOATIIET  I*2.oD-5 

TIME  «  TImE+El 

^F*’^KOUN‘f^io.  NUMBER  I  GO  TO  300 


XO  .■  X 
VO  «  Y 


UC  ■  U 
VC  -  V 
AOOOLD 
GO  TO  200 

s  t 

Mmtiiitimimumitmimmmmnsmmsmtiimmm 


AOD 


10 

8 

40 

I! 

80 

90 

100 

110 

120 

130 

140 

ISO 

160 

il°0 

190 

200 

82 
230 
240 
250 
260 
270 
280 
290 
300 
210 
320 
33  0 
’40 
-SO 
360 
370 
380 
390 
400 
410 
420 
430 
440 
450 
460 
470 
480 
490 
500 
510 
520 
530 
540 
550 
560 
570 
560 
590 
600 
610 
62  0 
630 
640 
650 
66  0 
670 
680 
690 
700 
710 
720 
730 
740 
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300  Sri riN< f, 500 j  nn.time 

WRITE  16.600)  AA, SUM, SUMA.X.V, LOSC, MMMNM 
600  CONTINUE 


ITERATIONS  IS 


500  FCAMAT  (IX, 'TOTAL  ELAPSEC  TIME  AFTER 
600  FORMAT  ( IX, IP  5c 20 .5,2 1  0 , /) 

?M> 

LBROUTINE  CURVE  (A.P.X.V) 

MPLICIT  R  E  AL*6  (A-h,0-2) 

X  •  A-(P*PI 
»  ■  P 
RETURN 
EKO 

$18 ROUT  INE  INCREP  (P.OELP) 

IMPLICIT  R  EAL*d  (A-H.0-2) 

P  ■  P*OELP 

RETURN 

ENO 
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SECTION  A-2 :  THE  STEEPEST  DESCENT  CONTOUR 

USER  INSTRUCTIONS 

This  is  the  FORTRAN  IV  program  for  the  numerical  inver¬ 
sion  of  the  Laplace  transformation  of  a  function  f(s),  which 
performs  contour  integration  along  the  steepest  descent  path 
in  the  complex  plane.  The  following  instructions  are  pro¬ 
vided  to  assist  the  user  in  familiarizing  himself  with  the 
program  so  that  he  can  adapt  it  to  his  particular  require¬ 
ments  . 

1.  SUBROUTINE  XMARCH  is  a  subroutine  which  tracks  a  speci¬ 
fied  v=0  contour  (steepest  descent)  in  the  complex  plane, 
once  this  contour  has  been  intersected  for  the  first  time. 

As  shown  in  the  program  listing,  the  following  two  equations 

X  =  XHOLD  -  l.D-l*DCOS( THETA) 

Y  =  YHOLD  +  l.D-l*DS IN (THETA) 

march  the  numerical  integrations  in  the  favorable  direction 
along  the  contour.  The  equations  are  parametric  in  THETA, 
which  is  the  search  angle  provided  by  SUBROUTINE  ROTATE. 

The  l.D-1  is  a  polar  radius  which  may  be  altered  to  a  larger 
or  smaller  increment  by  the  user,  or  accepted  as  it  appears. 

2.  SUBROUTINE  ROTATE  increments  the  angular  search  angle 
THETA  in  order  to  locate  a  new  point  on  the  v=0  contour.  The 
angle  is  presently  incremented  by  5.D-2  radians  with  each 
call  to  the  subroutine.  This  may  be  altered  as  required. 

3.  SUBROUTINE  VALUE  is  a  user  supplied  subroutine  which  cal¬ 
culates  the  real  and  imaginary  components  for  g(s)  =  exp(st) 

x  f(s)  ■  u(s)  +  iv(s).  For  examples  of  preparation  of  this 
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subroutine,  the  user  is  referred  to  Appendix  B. 

4.  The  mandatory  input  variables  to  the  main  program  are 
listed  as  follows: 

(a)  AA  -  REAL*8  starting  position  (s=AA)  on  the  real  axis 

(b)  AL  -  REAL*8  value  of  any  constant  associated  with 

the  function  f(s).  (Note:  If  the  function  f(s) 
involves  more  than  one  constant,  the  calling 
statements  to  SUBROUTINE  VALUE  must  be  modified 
accordingly. ) 

(c)  DELY  -  REAL*8  increment  of  the  point  spacing  along 

the  Bromwich  contour.  (Note:  to  alter  the  point 
spacing  along  the  v=0  contour,  the  user  is 
referred  to  1  above.) 

(d)  NREG  -  INTEGER*4  input  value  which  specifies  the 

first,  second,  third,  etc.  v=0  contour  above 
the  positive  real  axis,  which  is  to  be  followed 
in  the  numerical  integration  process. 

(e)  T  -  REAL*8  time  for  which  F(t)  is  to  be  evaluated. 

5.  Input  parameters  which  may  be  altered  by  the  user,  at 
his  discretion,  are  listed  below  as  follows: 

(a)  EPS  -  REAL*8  tolerance  parameter  for  the  termina¬ 

tion  of  numerical  integration.  The  program 
listing  presently  contains  a  value  of  l.D-11. 

(b)  NUMBER  -  INTEGER*4  number  of  successive  addend 

contributions  less  than  the  specified  value  of 
EPS  for  which  numerical  integration  is  terminated. 
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6.  The  output  variables  which  are  received  from  the  program 
are  defined  below: 

(a)  MMMMM  -  the  number  of  calls  made  to  SUBROUTINE  VALUE 

(b)  SUM  -the  total  addend  contribution  for  the  numeri¬ 

cal  integration 

(c)  SUMA  -the  absolute  value  of  the  total  integrand 

contribution  for  the  numerical  integration 

(d)  EL  -the  computation  time  in  seconds  required  for 

the  process 

7.  A  sample  program  output  is  as  follows 

694  MMMMM 

0. 724724D+0  0.324972D+1  SUM, SUMA 

0 . 707107D+0  0 . 326734D+1  SUM, SUMA 

4027  MMMMM 

0. 71090D+1  EL 


at  end  of  vertical 
section 
at  end  of 
computation 
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GLOSSARY  OF  VARIABLE  NAMES 

VARIABLE  DESCRIPTION  LEGEND 

AA  starting  position  on  the  real  axis  of  2 

the  Bromwich  contour  segment  of  the 
integration 

ADD  computed  value  of  the  addend  of  numeri-  1 

cal  integration 

ADDA  absolute  value  of  the  addend  of  numeri-  1 

cal  integration 

AL  constant  associated  with  f(s)  2 

CHECK  flag  to  ascertain  when  a  v=0  contour  1 

crossing  has  occurred 

DELY  increment  of  point  spacing  along  the  2 

Bromwich  contour  segment  of  numerical 
integration 

EPS  tolerance  parameter  for  termination  of  4 

numerical  integration 


EXAM  flag  to  ascertain  when  a  contour  crossing  1 

has  occurred 

J  integer  counter  used  internally  for  test-  1 

ing  contour  regions 


KOUNT  integer  counter  utilized  to  terminate  1 

the  numerical  integration  process 

KSIGN  integer  counter  utilized  to  record  the  1 

number  o  f  v=0  contour  levels  crossed 

LEVEL  integer  flag  assigned  to  regions  of  1 

positive  and  negative  v 

MMMMM  integer  counter  for  the  number  of  calls  3 

made  to  SUBROUTINE  VALUE  during  each 
segment  of  the  numerical  integration 
process 

NREG  the  desired  v=0  contour  level  above  the  2 

positive  real  axis  which  is  to  be  tracked 

NUMBER  integer  input  for  the  number  of  succes-  4 

sive  addend  contributions  less  than  the 
tolerance  parameter  required  for  termi¬ 
nation 
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SUM  total  addend  contribution  for  the  3 

numerical  integration  process 

SUMA  absolute  value  of  the  total  addend  3 

contribution  over  the  contour  of 
integration 

T  time  for  which  F(t)  is  to  be  evaluated  2 

TERM  termination  parameter  for  the  absolute  1 

value  of  the  addend 

TEST  Flag  to  ascertain  when  a  contour  1 

crossing  has  occurred 

THETA  search  angle  used  in  the  contour  1 

tracking  algorithm 

THETAO  previous  value  of  THETA  1 

THNEW  predicted  angular  location  of  a  new  1 

point  on  the  contour  from  a  known 
previous  point • 

THONE  coordinate  utilized  in  linear  inter-  1 

polation 

THTWO  coordinate  utilized  in  linear  inter-  1 

polation 

U  real  part  of  g(s)  1 

UHOLD  previous  value  of  U  1 

V  imaginary  part  of  g(s)  1 

VHOLD  previous  value  of  V  1 

VNEW  actual  value  of  the  imaginary  part  of  1 

g(s)  for  the  predicted  v=0  location 

VOLD  previous  value  of  V  1 

X  real  part  of  s  1 

XHOLD  previous  value  of  X  1 

Y  imaginary  part  of  s  1 

YNEW  predicted  elevation  of  the  desired  v=0  1 

contour  level 
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LEGEND 


1.  No  action  required  by  the  user 

2.  Mandatory  user  input 

3.  Appears  as  program  output 

4.  May  be  altered  by  the  user  at  his  discretion 
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CALL 

CURVE  AND  VALUE 


f 

COMPU' 
VALUE 
INTEGRi 
ITS  AE 
VAJ 

— 

rE  THE 

OF  THE 

AND  AND 
ISOLUTE 

LUE 

_ y ... 

SUM  THE 
CONTRIBUTIONS 
OF  THE 
INTEGRAND 
AND  ITS 
ABSOLUTE  VALUE 


COMPUTE  THE 
PRODUCT  OF  THE 
OLD  AND  NEW  V 


73 


75 


LT. EPS 


IF 

DABS(V) 


GT.EPS 
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COMPUTE 
TERMED ABS ( ADD ) 


THIS  PROGRAM  COMPUTES  THE  INVERSE  LAPLACE  TRANSFGRM  OF  A 
FUNCTION  USING  THt  STEEPEST  DESCENT  CONTOUR 
IMPLICIT  RE AL*S  (A-H.O-Z) 

CALL  SETIME 

CALL  ERRSET  (207,256,-1,1.1,209) 

****#***********»%*****♦**********♦*****♦************************ 

* 

*  THIS  SEGMENT  OF  THE  PROGRAM  NUMERICALLY  INTEGRATES  ALONG  THE 

•  BROMWICH  CCNTOUR  UNTIL  IT  INTERSECTS  WiTrl  THE  DESIRED  V*0 
•CONTOUR.  LINEAR  I NTERPCLATI  CN  IS  USEO  TO  REFINE  TFE  LOCATION. 

TWOPI  •  4.U*0*OARSIN( l.C*0)  LE'oLNU 

T  »  TWOPI - M-  2 

KKKKK  «  0  ^ 

AL  »  1  .250—1- — - - _____ - 2 

J  «  1 
MMMMM  -  0 
SLM  ■  0.0*0 

sum  -  o.o*o 

AA  *  3.C-1 - _ 

X  -  AA 

V  »  0.0*0 

Y  FOLD  *  0.0*0 

OELY  *  2.50-3 - ^ 

NREG  *  3 _ ____ _ _ _ >- 

KS1GN  *  0 -  - 

LEVEL  •  -l 

CALL  VALUE  (X,Y,T,AL,LHCLD,VH0LD, MMMMM) 

200  Y  *  Y *OELY 

CALL  VALUE  (X, Y , T ,AL, L , V .MMMMM) 

J  ■  J  *1 

IF  (J.EC.2)  SIGN  •  V/CAESIV) 

IF  (J.GT.2)  TEST  *  V*SIGN 
IF  (TEST  .LE.0  .0*0  )  30  TO  400 
300  ACO  -  ( (U*UHOLJ)*(Y-VhCLC))/TWOPI 
SLM  •  S  UM* AOO 
ACOA  -  CABS (AOO  I 
SUMA  -  SUM  A*  AOO  A 

IF  (KKKKK. EC. 1)  WRITE  (6,1800)  AOO.ADOA 

UHOLO  •  U 

VHOLD  -  V 

YHOLD  ■  Y 

GC  TO  200 

400  KSIGN  «  KSIGN*1 
LEVEL  -  -LEVE,. 

IF  (KSIGN. eC.NREG)  GO  TC  500 
SIGN  •  -SIGN 
GC  TO  300 
500  YCNt  ■  Y-OELY 
YTVO  »  Y 
VCNE  •  VHOLO 
VTWO  •  V 

600  VNEW  *  ( VONE*(YTWO-YONE)/(VONE-VTWO))*YCNE 
CALL  VALUE  (X ,Y ME W, T , AL , UNEW, VN EW , MMMMM  ) 

IF  (OAbS(VNtW).LT.l.fi-ll)  GO  TO  9C0 


IF  (OABS(VNtW).LT .1.0-1 
TEST  *  VNEW*SIGN 
IF  (TEST. GT  .0.0*0)  GO  T 
VTWO  •  YNEW 


-11)  GO  TO  9C0 
TO  700 


VTWO  ■  VNEW 
GO  TO  800 
700  YQNE  •  YNEW 
VONE  ■  VNEW 
800  CONTINUE 
GO  TO  6C0 

900  ACO*  { (UN EW*UHOLO)*( YNEW- YHOLO)) /TWOPI 
ACOA  •  CABS ( AGO ) 

SLM  •  SLM* AOO 
SUMA  *  SUM A*AOOA 

IF  ( KKKKK.E C. 1 )  WRITE  (6,1800)  AOC.AOOA 
WRITE  (6,2100)  MMMMM 
WRITE  (6,2000)  SUM, SUMA 

* 


st$sss<ssssss«sssss$ss$s$s$ss$$$ss$$st$»s$$$s$s$ss$$*s$$ss$$$$s$s 

$  THIS  SEGMENT  OF  THE  PROGRAM  PERFORMS  NUMERICAL  INTEGRATION 
(ALONG  The  V  «  0  CONTOUR.  LF'3tl[\lC 

KOUNT  -O 

NUMBER  »  5 - -  <4 


EPS  -  I.  0-1 1 - - - - - — 4 

T FETA  »  C.DhO  ^ 

Y  -  YNEW 
J  -  I 
XFOL 0  »  X 
YHOLO  •  Y 
UHOLD  »  UNEW 
VHOLJ  »  VNEW 
MMMMM  *  0 
J  *  J+l 

CALL  XM  ARC  F  (XHOlO. YHOLO, THETA,  X,  Y) 

CALL  VALUE  IX  .Y , T , AL, U , V , MMMMM ) 

IF  IJ.E0.2I  SIGN  «  V /DAB  SI  V J 
IF  (J  .GT.2  J  EXAM  *  V  *  S  I  £N 
IF  (EXAM. LE.O. 0+0)  GO  TC  1200 
VOLO  *  V 
THETAO  «  THETA 

CALL  ROTATE  (LEVEL, SIGN, ThETA) 

GO  TO  1IC0 
ThCNE  «  THETAO 
THTWO  *  THETA 
VO NE  ■  VOLO 
VTWO  «  V 

THNE  W  »  I  I  (THON£-THTWC)/(VTWO-VCNE) )*VONE I +THONE 
CALL  XMARCH  I Xri 3L D, YHOLO , THNE W,  X , Y ) 

CALL  VALUE  (X ,Y ,T , AL , UN EW,VNEW , MMMMM J 

IF  (OABSI  VNEMJ.lT. 1.0-1U  GO  tj  itOO 

CHECK  *  VNEW+SIGN 

IF  (CHECK. GT. 0.0+01  GO  TO  1400 

THTwO  *  THNE W 

VTWO  *  VNEW 

GO  TO  1S00 

THONE  ■  THNEM 

VCNE  »  VNEW 

CCNTI  NUE 

GO  TO  1300 

ACO  *  (  (U.4EW+UHJLGi*<Y-YF.OLO)+<VNEW+VHOLO)*<X-XHOLD))/TWOPI 
ACOA  ■  0A8S( AUDI 
SUM  «  SUM+AOO 
SUMA  »  SUM A+ADOA 

IF  (KKkKK.EO.1)  WRITE  <6,18O0>  ADC, ADDA 


THIS  IS  THE  TERMINATION  CRITERION  SEGMENT  OF  THE  ALGORITHM. 


0000000000000*04*0 044004 *««***« 

0  THIS  IS  THE  TERMINATION  CRIT 

• 

TERM  *  DABS (ADO) 

IF  (T  ERM.LT. EPS  J  KOUNT*KOUNT +1 
IF  ( TERM. UE .EPS)  K0UNT*0 
IF  (KOUNT. EQ. NUMBER)  CO  TO  1700 

• 


THETA  »  THNEW 
EXAM  ■  -EXAM 
GO  TO  1000 

1700  WRITE  (6,2000)  SUM, SUMA 
WRITE  (6,2100)  MMMMM 
CALL  GET (ME  < I ET ) 

EL  ■  OFLOAT(IET) *2.60-5 
WRITE  (6,  1900)  EL 

»  S 

mmmsiiitsmtitiiimtisittmiiimmmuimttmmim 


*1900  FORMAT  Iu^IlaPSE^TIME  OURING  THE  COMPUTATION  IS 


17  UW  ■  VJ  '  ^  '  I  •**  T  ■  — 

2000  FCRMAT  <lX.2E20.oi 
2100  FORMAT  (  IX. 18) 

fuBROUTINE  XMARCH  (XHCLC .YHOLO , THET A ,X , Y 1 

IMPLICIT  REAL*8  (A-H.O-Z)  _ 

X  *  XHOLD-2 .50-3* CCOS (T H ETA h — — 

»  -  YHOLD+2. 50-3*0  SIN  (THETA)  — - 

RETURN 
ENO 


•  ,E20.5 J 

LtG 


E.N1 


100 


SUBROUTINE  ROTATE  (LEVEL .SIGN, THETA ) 

Implicit  real*8  ia-h,q-z) 

If  <  sf  GN-Llf  oio-*-q"?  TH|TAsTHETA-5.C-2 

H 


( SI GN.GT *0*0+0 )  THETA-THETA* 5.0-2 
I2i28?LT.O.O*OI  IHfTA»THETA*5.C-2 


IF  (SIGN.GT  .0.0*0)  Trt ET A*T HET A- 5 .0-2 
200  RETURN 
ENO 

LEGEND 


}490 
500 
510 
1520 
1530 
1540 
10 
20 
30 
40 
50 
60 

io° 

30 

40 

50 

$0° 

eo 

90 

100 


2.  Mandatory  user  output. 

4.  May  be  altered  by  the  user  at  his  discretion. 


5.  Polar  radius  R  =  .025  built  into  program,  may 
be  changed  by  user. 


SECTION  A- 3.  SHANKS'  ACCELERATOR 
USER  INSTRUCTIONS 

The  algorithm  of  this  section  is  very  similar  to  that 
of  Section  A-l.  Only  those  features  which  differ  have  been 
discussed  in  this  instruction. 

1.  Contour  integration  is  performed  using  a  simple  para¬ 
meterized  curve  for  a  finite  number  of  points  on  this  con¬ 
tour.  This  finite  number  of  terms  must  be  prescribed  by 
the  user.  The  integer  variable  M  is  used  for  this  purpose. 

2.  The  dimension  statements  at  the  beginning  of  the  program 
must  provide  storage  space  in  linear  arrays  which  are  equal 
to  or  larger  than  the  value  assigned  to  M.  (These  quanti¬ 
ties  will  subsequently  be  used  with  Shanks'  accelerator. 

Note  particularly  that  the  array  is  douoly  indexed.) 

3.  The  integer  variable  N  specifies  the  number  of  iterates 
of  Shanks'  accelerator  which  are  desired  by  the  user.  How¬ 
ever,  the  assigned  value  of  N  cannot  be  greater  than  6  as 
the  program  is  now  written. 

4.  A  sample  program  output  appears  in  Chapter  VI  of  this 
thesis.  Column  1  represents  the  natural  result  obtained 
after  M  evaluations  of  the  addend,  without  application  of 
the  accelerator. 


no  onoonooo  noooonm 


THIS  PROGRAM  CALCULATES  THE  INVERSE  LAPLACE  TRANSFORM  OF  A 
FUNCTION  USING  A  FINITE  NUMBER  OF  TERMS  ALONG  A  SIPPLE 
PARAMETERIZED  CONTOUR  CF  INTEGRATION  <UTH  SHANKS'  ACCELERATOR 
IMPLICIT  REAL *8  1A-H.O-Z) 

DIMENSION  AUO.IJO)  ,  ACOUOO) 

OIMENS1  CN  XX<  10J)  ,YYl 100)  ,UU( 100)  ,VV(10  0) 

CALL  ERRSET  I  20 7, 256,- I.  1. 1. 20S) 

MPMMM  *  0 

TWOPI  »  4.0+0*0ARSIN<1.0+0) 

AL  «  1.250-1 
T  »  TWO PI 
AA  •  3.0-1 

P  -  0.0+0  ,  —  " 

M  «  20  __ 

NN  •  0 
N  «  5 

CALL  CURVE  <AA, P.XO.YO) 

CALL  VALUE  I  XO, Y0 , T , AL , IQ, VO , MMMMM) 

t»$S$$$SS$S$»S$$S$S$S*S$S$$$SSS«S*SS$SSSSS$$$$*S$SSS$SSS$S$S$SS$SS 

s  * 

»  THE  FCLLGWI.NG  SEGMENT  OF  THE  PROGRAM  PERFORMS  THE  CONTOUR  * 

AINTEGR ATION  FOR  A  FINITE  NUMBER  CF  TERMS  USING  THE  SIMPLE  » 

*  PARAMETERIZE!)  CURVE.  i 


$  * 
200  NN  «  NN  +  1 

CALL  INC  PEP  (PI 
CALL  CURVE  (AA,P,X,Y) 

CALL  VALUE  (X, Y.T.AL.U.V, MMMMM) 

ACO(NN)  *  UV+VO)*(X-XO)  +  (U+JO)»(Y-YOn/TWOPI 

XO  «  X 

YO  -  Y 

UC  «  U 

VO  «  V 

XXI  NN)  -  X 

YY(NN)  -  Y 

UU(NN)  *  U 

VV(NN)  *  V 

IF  INN.EO.M  )  GO  TO  300 
GO  TO  200 
300  CONTINUE 

*♦♦******»*****♦♦*******♦*******♦********»**♦**•■?*****♦♦♦♦♦♦♦♦**»«« 
*  * 

♦  THIS  SEGMENT  OF  THE  PROGRAM  IMPLEMENTS  SHANKS'  ACCELERATOR  IN  * 

♦  ACCORDANCE  WITH  EQUATION  28  OF  CHAPTER  6  OF  THIS  THESIS.  * 

♦  * 
8  *  1.0+0 

EPS  *  1.0-12 


OC  400  1-1, N 


OC  400  J-l ,  M 
400  A 1 1 , J )  -  0.0+0 


00  500  J-t ,M 
JM  -•  J-l 
B  -  AOO(J) 

IF  IJ.EC.l)  A 1 1  , 1 )  -  e 

IF  (J.GT.l)  A  (  1 ,  J  )  -  B+AU.JM) 

CCNT1NU6 
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?8  S°?-i-2’N 

K  •  2*1-1 


00  600  J«K,M 
JP  -  J-l 
JHM  -  J-2 
AJ l.JI  -  1.0*0 

DEN  -  A(IH,  J)+A(IM,JMH)-2.C+0*A(IP,JM) 

It  <OAB$(DEN)  .ST.  EPS)  All  ,  J)  -  ( A(  I  P ,  J)  *A  < 1 P ,  JMP)  -A<  I  P.  JM  J**2  J /OEN 
600  CONTINUE 


OC  700  J*2 >  P 

700  WRITE  (6.800  J  .  (  A(  I .  J  )  .  1*1  ,N  J 

*  * 
****************************************************************** 


800  FORMAT  (5X, I3.6E18.14 J 
END 

SUBROUTINE  CURVE  (A.P.X.Y) 
1PPUCIT  RE  AL*8  4A-H.O-2J 
X  «  A-4P*P» 

Y  -  P 

RETURN 

ENO 

SUBROUTINE  INCREP  (P) 

IPPUCIT  REAL*8  4A-H.G-21 

P  -  P  +  1.0-1 

RETURN 

ENO 
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SECTION  A-4.  SPECIAL  CONTOUR  ALGORITHM  FOR  CASES  AS 
TREATED  IN  CHAPTER  EIGHT 

USER  INSTRUCTIONS 

The  integration  contour  used  in  this  section  is  a  vertical 
path  to  a  specified  value  of  y,  follov/ed  by  a  simple  parabola. 
The  algorithm  of  this  section  is  very  similar  to  that  of 
Section  A-l.  Only  those  items  which  differ  are  discussed  in 
this  section. 

1.  Contour  integration  is  performed  in  the  initial  segment 
of  the  program  along  the  Bromwich  contour.  The  integer  vari¬ 
able  JJJJJ  specifies  the  number  of  poles  of  the  function 

1  r<c 

f (s)  =  —  tanh(-=-)  above  the  real  axis  that  are  to  be  enclosed 

S  A 

to  the  left  of  the  distorted  Bromwich  contour,  which  is  used 
in  the  second  segment  of  the  program. 

2.  The  variable  HALT  computes  the  elevation  above  the  real 
axis  corresponding  to  the  input  value  of  JJJJJ  for  the 
specified  hyperbolic  tangent. 

3.  The  variable  FLAG  is  internally  assigned  a  value  of  zero 
or  one,  corresponding  to  the  first  or  second  segments  of  the 
contour  of  integration  which  the  algorithm  employs.  No 
action  is  required  by  the  user. 

4.  The  integer  variable  KSTART  is  utilized  in  the  outer  loop 
of  the  program  to  increment  the  value  of  t  for  which  F(t) 

is  computed. 

5.  The  linear  arrays  STORE  and  STASH  are  used  to  convert  the 
values  of  f(t)  and  t  to  single  precision,  which  is  required 
for  the  graphics  package  associated  with  the  IBM  360/67. 


6.  SUBROUTINE  VALUE  contains  the  real  arithmetic  coding  of 
the  function  g(s)  for  the  hyperbolic  tangent  function  of 
Chapter  VIII  of  this  thesis.  Standard  trigonometric  identi¬ 
ties  for  the  hyperbolic  functions  of  a  complex  argument  have 
been  utilized. 

Note:  This  procedure  is  obviously  keyed  to  the  particular 
function  f(s)  =  ^  tanh(^).  The  user  must  modify  it  for 
other  functions.  One  obvious  modification  which  would 
generalize  the  procedure  would  simply  be  to  specify  rather 
than  calculate  the  value  of  HALT,  i.e.,  input  HALT  rather 
than  JJJJJ . 
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THIS  PROGRAM  COMPUTES  THE  INVERSE  LAPLACE  TRANSFCRM  OF  A 
FUNCTION  US  IN G  THE  SIMPLE  PAR AM£ TERI ZEO  cONTOUR  OF  INTEGRATION 
WITH  MOOIF ICATI ON  TO  VAST  THE  tLEVATION  FROM  WHICH  THE  CURVE  IS 
LAUNCHED.  THE  BROMWICH  CONTOUR  IS  USED  IN  THE  INITIAL  SEGMENT. 
IMPLICIT  REAL *8  (A-H.O-Z) 

REAL*4  STORE,  STASH 
DIMENSION  STORE! SO),  STASHI 50) 

CALL  ERRSET  1/07,23 6, -1,1,1,209) 

TWOPI  -  4.  D40*0AR SI N( 1.0*0) 

AL  •  1.041 
JJJJJ  «  15 
KKK«  -  0 
NUMBER  -  5 
KST ART  «  50 


DO  400  KKK  » 
FLAG  *  0.040 


1,  K START 


HALT  *  IDFLOATI  JJJJJ) *TW0PI)/(2.D40*AL) 

ACD  »  0.040 

ACOA  -  0.040 

SLM  «  0.040 

SUMA  -  0.040 

LCSC  ■  0 

MMMMM  ■  0 

AK  «  KKK 

AA  a  3.0-1 

T  •  5. 0- 14AK* 1.04 0 

P  a  0.040 

KCUNT  ■  0 

NN  a  0 

EPS  -  i.O-ll 

OELP  ■  1.0-1 

CALL  CURVE  ( AA, P, HALT ,FLAG ,X0 , VO) 

CALL  VALUE  (XO,YO  ,T  ,AL,UC,  VO,  MMMMM) 

TIME  ■  C.D40 
CALL  SETIME 
200  NK  ■  NN4 1 

CALL  /NCR  EP  <  HALT,  0£  LP.F.FLA5  ) 

CALL  CUFVE  (AA, P, HALT , FLAG, X, V ) 

CALL  VALUE  (X,Y,T,AL,U,V, MMMMM) 

AOO  •  ((V4V0)*(X-XJ)4(U4U0>*(  Y-VO))/TWOPI 

IF  (NN.EC.l)  AOJOLOa AOC 

PROD  -  AOO*  AOOOLJ 

IF  (PAOO.LT. 0.J40)  LO  SC-LOSC* 1 

AOOA  -  OABS (ADO) 

SUM  >  SUM4A00 
SUMA  >  SUM  A  4A0DA 
TEST  »  OABS I AOO) 

IF  (TEST.LT.EPS)  K0UNT»KCUNT4 1 
IF  (TEST.GE.EPS)  KOUNT «0 
CALL  GETIME  (IET) 

EL  a  OF  LO AT ( I ET ) *2.60- 5 
TIME  »  TIME*EL 
CALL  SETIME 

I F(KOUNT.EQ. NUMBER)  GO  TO  300 
XC  a  x 
TO  -  Y 
UO  a  U 
VC  a  V 

AOOOLO  a  AOO 
GO  TO  200 
300  CONTINUE 

WRITE  (6,500)  NN.TIME 

WRITE  (6,600)  T,  SUM,  SUMA ,  X,  Y,  LOSC,  MMMMM 
STORE(KKK)  *  SNGL(T) 

....  STASH  (KKK)  »  SiNGL  (SUM) 

400  CONTINUE 

CALL  PL CTG (STORE, ST ASF, SO, 1,1,0, • T • , 1,  »F ( T) ', 4 ,0.0, 55. 0,-2. 0, 2. C, 
*!ift6PL&T (0.0,0.0,949) 


500  FORMAT  (IX,  ‘TOTAL  ELAPSEO  TIME  AFTER  M5,*  ITERATIONS  IS  *,E20.5) 


600  FORMAT  ( XX ,  1P5E20  .5, 2  1 8,  / ) 


llBROUTl 


100  X 
2  00 


.  INS  CURVE  IAA, F, HALT, FLAG, X.Y) 
PPLICIT  RSAL*a  (  A-h,  0-Z  ) 

F  IP. LT. HALT. ANU. FLAG. AS. 1.0+0)  GO  TO  100 
■  A A- 1  IP- HALT )* IP- HALT ) I 

GC*TO  200 
AA 

Y  ■  P 
RETURN 
END 

SUBROUTINE  INCREP  I  HALT, JELP, P ,FL AG ) 
IMPLICIT  R6AL*3  I A-H, 2-21 
P  »  P*OELP 

IF  IP. EQ. HALT)  FL AG*  1.3  +  0 

RETURN 

END 

SUBROUTINE  VALUE  I  X, Y , T  ,  AL, U, V, MMHMM) 
IMPLICIT  RE  AL*3  1A-H.C-2) 

MMMMM  «  MMMMM+i 
I  AL  *X  ) 

I  AL*Y ) 

OSINHI TEML) 

CS I N IT EM2) 

0CCSH(TcMl)+3CCSITEM2) 

TEM3/TEM5 
TEM4/T  EM5 
T  I  ■  I  TEM6*X)  +  1  TEM7*YJ 
T  *  ITEM7*X)-I  TEM6*Y) 

T  0  *  IX*X)*I7*Y) 

T  1  •  TEM8/TEM10 
~  _I  TEM9/TEM10 
TEM13  -  PEXPI X*T)*OCC«I Y*T) 

TEM14  »  L6XPI X*T) *JSINI  Y*T) 

U  -  !TEM13*TEMU)-<TEMU*TEM12) 

V  *  ITEM14*TEMU)*ITEP13*TEM12) 

RETURN 

ERO 
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OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 


APPENDIX  B 


LISTING  OF  TEST  CASES  AND  EXPERIMENTAL  RESULTS 

The  Laplace  transformation  pairs  f(s)  and  F(t)  which  we 
have  used  as  test  cases  are  listed  within  this  appendix. 

Also  included  for  each  function  f(s)  is  the  listing  of  the 
SUBROUTINE  VALUE  and  a  typical  result  of  numerical  inversion 
at  one  specific  value  of  t.  Case  9  also  shows  a  plot  of 
F(t)  versus  t  obtained  numerically. 

The  results  obtained  by  numerical  inversion  have  been 
compared  where  applicable,  with  tabulated  values  listed  in 
(1).  In  all  other  cases  the  value  of  F(t)  has  been  cal¬ 
culated  on  the  IBM  360/67. 

The  symbols  which  are  used  in  the  presentation  of 
results  are  defined  below: 

AL  -  constant  associated  with  f(s) 

BE  -  constant  associated  with  f(s) 

T  -  time  for  which  F(t)  is  evaluated 

AA  -  real  axis  starting  position  of  the  integration  contour 

xf  -  final  value  of  x  at  termination 

yf  -  final  value  of  y  at  termination 

MMMMM  -  total  number  of  calls  made  to  SUBROUTINE  VALUE 

LOSC  -  number  of  sign  changes  between  successive  addends 

SUM  -  value  of  F(t)  (equation  (16)) 

SUMA  -  value  of  G(t)  (equation  (19)) 

TIME  -  clock  time  (sec)  required  for  numerical  integration 
ANALYTICAL  -  analytical  value  of  F(t) 

In  all  cases  we  used  N  =  5  and  EPS  =  l.D-11. 
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l/(s  -aJ)  51  /s-a 


Case  1 


c 


IBStfJR  SlSJlxiii  tft 


U,v  .MMvmm) 

01  OF  LAPLACE 


TRANSFORMS  8Y  SPIEGEL 


MMMMM 
X  PTR  ■ 
XPTI  - 
CXPT  » 
CMJM  » 
r.OENOM 
CRESLT 

u 

V 

RET'JRN 

EIO 


UMVIMM  +■  1 

X*T 
Y*T 

OCM’LXIXPTR.XPTI  ) 
COE  XPICXPT) 

*  DCM?LX(A,Y) 
t,  *  CMUM/COEMOM 
OREAUCRESLTj 

0 1  MAGI  C  RE  aLT) 


AL  «  1.25D-1 
T  »  8.D+0 
AA  «  0.2 
Xf  -  -U.6U. 

y F  *  2.2 
MMMMM  *  23 
LOSC  ■  5 

SUM  -  l.OOOOOD+O 
SUMA  ■  1.42588D+0 
TIME  «  0.7979^0-1 
ANALYTICAL  “  l.OOOOOD+O 
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Case  2 


SUBROUTINE  VALUE  ( X» Y, T • AL » U, V , MMPMM) 

THIS  IS  "A8LE  EVTRr  NUMBER  02  OF  LAPLACE  TRANSFORMS  8Y  SPIEGEL 
IMPLICIT  REAL*8  < A,3,0-H,0-Z) 

IMPLICIT  COMPLEX* 16  Id 
MMMMM  «  MMMMW  «■  1 
CS1  *  DCMPLX(X.Y) 

XPTR  =  X*T 
XPTI  *  Y*T 

CX»T  *  OCMPLXI  XPTR, XPTI ) 

CNUM  =  COEX  R( CX  3T  ) 

COENOM  *  OCMPLX (X ,Y> 

COENOM  »  C  S1*CS1 
CRESLT  =  CN'JM/ COENOM 
U  *  ORE  ALICSESLTJ 
V  *  DIM AG(CRESLT) 

RETURN 

EM) 


AL  =  1.25D-1 
T  *  3.D+0 
AA  -  0.2 
XF  »  -4.64 

yF  *  2-2 

MMMMM  *  23 
LOSC  ■  4 

SUM  *  8 . OOOOOD+Q 
SUMA  «*  8.04195D+0 
TIME  -  0.53248D-1 
ANALYTICAL  *  8.00000D+0 
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Case  3 


S  UBROUT  INE  V  ALU  6  (X ,  Y  ,  T  ,  it , N^W S=  ,  U.  V,  MMMMM  ) 

THIS  IS  TABLE  ENTRY  NUMBER  03  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMOLICIT  °S1L*8  (  A. 3  ,0-H,0-Z  ) 

IMPLICIT  COMPLEXES  ICi 


MMMMM  •  MMMMM  ♦  1 
ZERO  *  0.0+0 
CS1  «  OCMt>LX(X,YJ 
XPTR  »  X*T 
XPTI  »  Y*T 

CXPT  *  DCMPLXIXPTR.XPTI) 
CHUM  ■  CDE XP( CX°T ) 

CDENOM  »  CSIM*N»OWER 
CRESLT  -  CNUM/CDENOM 
U  •  DR  E  &L  (  CRESL  T) 

V  -  DIM AGICRESLT I 

RETURN 

ENO 


AL  -  1.25D-1 
N POWER  -  3 
T  «  8.D+0 
AA  •  0.3 
Xy  ■  -4.54 
yp  -  2.2 
MMMMM  »  23 
LOSC  -  3 

SUM  ■  3.20000D+1 
SUMA  -  3.20008D+1 
TIME  «  0. 8369^-1 
ANALYTICAL  *  3.20000D+1 
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Case  4 


Wy’K'W.smrH’HShl*  egei 

I WPLIC It  real *8  ( A, 3, C-h,l-Z ) 

IMPLICIT  C0MPLEX*16  (C) 

OIMEMSICM  NC  ALL  15) 

MMMMM  ■  MMMMM  ♦  I 

CS1  *  OCMPLXIX.YI 
XPTP  «  X*T 
X  PTI  ■  Y*T 

CXPT  •  DCMPLXUPTR.XPTI! 

CNUM  «  COEXPt  CXP7 ) 


CALL  CPCHEP  ICS1.COENCM, POWER, l.NCALL) 
CRESLT  *  CMUM/COENOM 

v :  8fSSb{8^I!tTTl 


retupn 

END 


AL  -  1.25D-1 
POWER  -  1.7D+0 
T  -  8.D+0 
AA  -  0.2 
XF  -  -4.64 
yF  “  2.2 
MMMMM  ■  23 
LOSC  ■  4 

SUM  ■  U.71815IHO 
SUMA  »  U.8324D+0 
TIME  *  0.85098D-1 
ANALYTICAL  »  4.71815D+0 
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Case  5 


C 


••PLICIT  P. E AC *8  (  A, 3,  C-l-»0-Z) 
MPtniT  COMPLEX*  16  (C) 

•  *****  *  1 
:R0  »  O.D+O 
II  »  DCMPLXIX.Y) 

TP  *  X*T 

ISt :  fe»w?:x'Tn 

£rESLt  «  CNUM/COENOM 


i 

M 

L 

XP 


PREAC|C 


U  _ - 

V  »  DIMAG 
RETURN 
END 


■SLT) 

•SLTI 


TRANSFORMS  BY  SPIEGEL 


AL  -  1.250-1 
T  ■  8.D+0 
AA  •  0.3 
XP  -  -4.99 
yP  •  2.3 
MMMMM  •  23 
LOSC  •  5 

SUM  -  2.71828D+0 
SUMA  •  3.59474CH-0 
TIME  •  0.890500+0 
ANALYTICAL  •  2.7182*D+0 


|gs»My 


Case  6. 


Ts0“s1*aU,LE^5},N;45U‘,oS“",ULl»L;CE"?»«F3»«  6Y  S.I EGEL 


••PL  1C  IT  PEAL  *8  M,  9,0-t-»CW) 
IMPLICIT  COMPLEX* 16  (C) 

HVIMMM  ■  MMMMM  ♦  l 

ZERO  ■  0.0+9 
C$1  »  OCMBLXl  XfY> 

CS2  »  OCMP  LX ( AL «  ZERO ) 

XPTP.  -  X*Z 
XPTI  ■  Y*T 

CXPT  •  OCMPLXIXPTR, XPTI ) 

CWH  •  CDEXPICXPTJ 
CS3  »  CSl  -  CS2 
CSENOM  »  C$3*»N,>OWER 
C RE SL T  ■  CNUM/COENOM 


OREA^CCRE  SLT) 


RETURN 
ENO 


AL  -  1.25D-1 
N POWER  »  3 
T  -  8.0+0 
•  AA  •  0.4 

XF  »  -4.44 
yF  ■  2.2 
MMMMM  *  23 
LOSC  ■  3 

SUM  ■  8.698500+1 
SUMA  -  8.69965D+1 
TIKE  ■  0.757900-1 
ANALYTICAL  »  8.69850D+1 
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Case  7. 


1  I 
«  1 


•»■]  =  *a 


SUBROUTINE  VALUE  IX  *Y  .  T  ,  AL.MPOWER  *U  «V  ,  MMM**M.NC  ALL  ) 

THIS  IS  TABLE  ENTRY  NU*3EP.  07  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 

iNPUCIT  REAL *8  (  A,9,C-F,0-Z) 

PPLICIT  C0MPLEX*16  (C) 

IMgNSION  VC  ALL (5) 


MMMMM  ■  MWMMM  *  I 
ZERO  -  0. 3+0 
C$1  ■  OCMPL  X( Xf  Y ) 

CS2  *  OCMPLX  ( ALtZERO) 

CS3  •  C  SI  -  CS2 
XPTR  -  X*T 
XPTI  ■  Y»'r 

CXPT  »  OCMPLX(XPTP.,XPTI) 
CN»M  -  COEXPICXOT  ) 


LTVin  »  cue  A  K  \  CA  **  I  | 

CALL  CPGWER  I CS3 . COENCM. POWE R,1 .NCALL) 
CPESLT  -  CNUM/CD6N3M 
'J  *  OBEALICRESLTI 
V  -  OIMAG(CRESLT) 

RETJRN 


AI>  ■  1.25D-1 
POWER  *  1.7D+0 
T  ■  8.D+0 
AA  ■  0.3 

Xp  .-4.54- 5* 

yT  • 

MMMMM  -  23 
LOSC  •  4 

SUM  ■  1.282530*1 
SUMA  -  1.305330*1 
TIME  ■  0.969280-1 
ANALYTICAL  -  1.282 530*1 
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CdSG  8 


SU8RTUT  If j H  VALJ6C  X,  Y»  Tf  Al,  U#  ) 

IdMitfrViiE.iTK.srafS-?? 0F  L1'L,CE  T"“sF,a,'s  «  s'iemi- 

I f ®L1 C IT  C1MPL  EX* 16  ( C  J 
HHMMM  «  ♦  l 

ZERD«0.0«-0 
CS»OC*PLX(X,YI 
CALOCNPLXt  SL .ZEftO) 

CT;OCMPLX(T,zeRO) 

COEN  -  <<C**CS)+<  CAl*CAU> 

CEXP«:OEXP(C<*CT) 

C  ■  C EX  P/COEN 
U-OPEAL(C) 

V-OINAGCC) 

RETURN 

END 


AL  ■  X.25D-1 
T  »  8.D+0 
AA  -  0.3 
Xp  -  -4.54 

yF  •  2.2 
MMMMM  -  23 
LOSC  -  4 

SUM  •  6.73177D+0 
SUMA  ■  7.03097D+0 
TIME  »  0.86294D-1 
ANALYTICAL  •  6.73177D+0 
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Case  9 


cot  at 


C 


SUBROUTINE  VALUE  (X,  V.  T,  AL,  U,  V,  MMM 
THIS  IS  TA8L E  ENTRY  NUMBE*  09  OF  LAPLACE 
IMPLItlT  PE AL*8  < A,b . D-H ,0-Z) 

implicit  complex-ho  «c ) 

MMMMM  -  MMMMM  ♦  I 
ZERO-O. D+O 
CS*OCMPLX(X,Y) 

;AL-0CMPLX(  AL.ZERO) 

:t«o:mplxit.zer3) 

;  OEN*C  S  **  2+C  AL  **2 

:<s/coen  _ 

:exp«cdexp<cs*ct> 

C»C*CEX  P 
U«OREAL(C) 

V*OIMAG(C) 

RETURN 
END 


HUM  Ml 

‘  TRANSFORMS 


BY  SPIEGEL 


T  ■  0.D+O 
AA  ■  0.3 
X?  -  -U.99 
yp  -  2.3 
MMMMM  ■  24 


LOSC  *  5 
SUM  *  5.40302D-1 
SUMA  •  2 . 00010D+0 
TIME  »  0.1113BD+0 
ANALYTICAL  •  5.40302D-1 
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Case  10 


IMPLICIT  R  £iL*8  < A,B.D-H,0-2» 
IMPLICIT  Cn^PLEX+lfr  (C) 


TRANSFORMS  BY  SPIEGEL 


IMPLICIT  Cn^PLEX*l&  (C) 

MMMMM  *  MMMMM  ♦  L 
Z  ER3  »  0.0*0 
CSl  «  DCMPLX(X,Y) 

CS2  •  ocMPLxcae.  zeroi 

CS3  *  DCMPLXI AL .ZERO > 

XPTR  ■  X*T 
XPTI  *  Y*T 

CXPT  ■  OCMPLX(XP*R,XPTI> 
CNUM^*r COE XP( CX PT ) 

COENOM*  »*( CS4+CS4 )+(CS3+CS3) 
C  RE$LT  *  CNiJM/COENOM 
U  -  OBEAUCRESLTJ 
V  -  01 “ AG (CPESLT ) 

RETURN 

END 


AL  »  1.25D-1 
BB  -  1.250-1 
T  ■  8.D+0 
AA  «  0.3 
Xy  «  -4.54 

yF  ■  2.2 

MMMMM  -  23 

LOSC  ■  4 

SUM  •  1.829880+1 

SUMA  •  1.88117D+1 

TIME  *  0.700440-1 

ANALYTICAL  ■  1.829880+1 
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K _ ^JL 


Case  1 1 


^r1 


»-> 

(•“  4)!  +  «» 


«**  cosat 


imi  SMiiW  H  ^Wll^'CSSi.  T«**».s  BY  SP I  £GE  i. 

lStlel?  MiKIxltiW’5*" 

MKMMM  m  MMMMM  ♦  1 

ZERO  *  0.0+0 

CSl  ■  DCM°LXU.YJ 

CS2  «  DCMPLXUE.ZERO) 

CS3  -  DCMPi.XUL.ZERO) 

X  PTR  »  X*T 

Ml  p25L-»ff  ItoSJi 'cxot I 

£BeNOMCJ1 ( CS4*CS4)+(CS0*CS3) 

CRESL"  *  CNtJM/CDE'MOM 
U  •  OREAL(CRESLT) 

V  »  DIMAOCRESLT) 

RETURN 

END 


AL  -  1.25D-1 

BE  -  1.25D-1 

T  -  8.D+0 

AA  *  0.4 

Xp  -  -4.89 

yp-  2.3 

MMMMM  «  24 

LOSC  «  5 

SUM  ■  1.46869D+0 

SUMA  ■  4.61766D+0 

TIME  ■  0.7S130D-1 

ANALYTICAL  •  1.46869D+0 
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Case  12 


SUBROUTINE  VALJ E( X , Y, Tf AL. U, V, ) 

THIS  IS  TABLE  ENTRY  MUMEER  12  OP  LAPLACE  TRANSFORMS  BY  SPIEGEL 
ILLICIT  R  £  AL  *  8  < A ,B , O-H.O-Z J 
IMPLICIT  C0«PLEX*16  ICi 


MMMMM  »  MMMMM  «■  1 
ZERT*O.D*0 
CS*OCMPLXIX,YJ 
CAL«OCMPLX(  AL  .ZERO) 
Ct»OCTPLX(T,ZERO) 

COEN  «  <  CCS*CS)-( CAL*CAL) ) 
CEXP«:OEXP(CS*C  T) 

C  ■  C EXP/COEN 

wnmw 

RETURN 

END 


AL  ■  1.25D-1 
T  »  8.D+0 
AA  «  0.3 
Xy  •  -4.5<F 

yF  -  2.2 

MMMMM  -  23 

LOSC  »  U 

SUM  ■  9.*»0l62D-*-0 

SUMA  ■  9.^**9DtO 

TIME  ■  0.10691D+0 

ANALYTICAL  ■  9.90162D+0 


Case  13 


coahal 


SUBROUT  IME  V ALJ EIX . Y . T ,  AL.U,  J 

THIS  IS  TABLE  ENTRY  &UM8ES  13  OF  LAPLACE 
IMPLICIT  RE4L*8  I  A* 3 .D-rt»3-Z ) 

IPPLI^IT  COMPL EX* lo  (Cl 


MMMMM 


MMMMM  ♦  i 


ZERO*  0.0*0 
CS»0CMPLX( X » Y 1 

:al»ocmplxi al.zero) 

CT«0CMP  LX ( T»  ZERO  1  „ 

, DEN  •  <CCS*CS1-(CAL*CAL)1 

>cs/:oen 

'  EXP*  CO  EX3  (  CS  *CT  ) 

_<*C*CEXP 
U«0RE  AL  ( C ) 

V«0IMAG(C» 

RETJRN 

END 


TRANSFORMS  BY  SPIEGEL 


AL  *  1.25D-1 
T  •  8.0+0 
AA  »  0.3 
Xp  «  -H.99 
yF  -  2.3 

MMMMM  ■  2d 
LOSC  -  5 

SUM  -  1.5*3080+0 
SUMA  •  2.52*1010+0 
TIME  -  0.85176D-1 
ANALYTICAL  »  1.5**308D+0 
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t „  _J»  - _ J-  JV 


Case  14. 


tilth  af 


SUBROUTINE  V AL0E< X,Y,T,AL,  BE,  U,V, MMMMM)  TDAuec_BM-  ftv  cot fgfl 
THIS  IS  TABLE  ENTRY  NUMBER  14  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  R  E  A  L  *  8  ( A,B ,0-rt,0-Z >  ”* 

IMPLICIT  CO*PLEX*U  (Cl 
MMMMM  ■  MMMMM  ♦  1 
ZEPO  *  O.D+O 
CS1  *  OCMPLX(X.Y) 

CS2  -  OCMPLX(BE,ZERO) 

CS3  «  DCM°LX< AL.ZERQ) 

XPTR  =  X*~ 

CXPT  *  CCMOLXIXPTR.XPTI  ) 

CNJM  «  COEXP(CXPT) 

CS4  »  CS1  -  CS2 

COENOM  «  (CS4*CS4)-(CS3*CrR) 

CRESf  *  CNtlM/COENOM 
U  *  ORE  AL  ( C3  ESlT) 

V  *  DIMAGICRESLT) 

RETURN 

ENO 


AL  -  1.25D-1 
BE  »  1.25D-1 
T  *  8.D+0 
AA  *  0.5 
Xp  -  -A. ?9 
yF  «  2.3 

MMMMM  «  24 
LOSC  *  4 

SUM  -  2.55562D+1 
SUMA  -  2.81818D+1 
TIMS  «  0.74932r-l 
ANALYTICAL  *  2.555620+1 


.  .I  i  -  I.  . 
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Case  15 


■*1  »->  — 

(t  -  *>=  -  a1  — 


**•  coih  at 


SUBROUTINE  VALUEU.Y.T.AL.RE.U.V.mmMMM)  , 

THIS  IS  TABLE  ENTRY  NUMBER  15  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REAL *8  ( A, 5 fO-HtO-Z ) 

IMPLICIT  COMPLEXES  (CJ 
NMMMM  *  MMMMM  ♦  1 
ZERO  -  O.D-t-O 
CSl  *  OCMPLX(X,Y) 

CS2  »  OCMPLXOE.ZERO) 

CS3  ■  0CMPLX(AL,ZEP.D) 

XPTR  ■  X*T 
XPTI  *  Y*T 

CXPT  »  OCMPLXIXOTR.XPTI) 

CNUM  -  ICSI-CS2)*C0EXP<CXPT) 

CSA  »  CSl  -  CS2 

tOENOM  -  ( CS4*C$*  >-( CS  3*C$3 ) 

CRESLT  «  CNUM/:06NOM 
U  ■  OP  EAL ( C°E  SL  T) 

V  ■  OIMAG(CRESLT) 

RETURN 

ENO 


AL  ■  1.25D-1 
BE  -  1.25D-1 
T  ■  8.D+0 
AA  •  0.5 
Xp  -  -4.79 

yF  -  2.3 

MMKMM  -  24 

LOSC  »  5 

SUM  »  4.19453D+0 

SUMA  •  9.63599D+0 

TIME  »  0.77740D-1 

ANALYTICAL  -  4.19453D+0 
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Case  16. 


-1 


(» - «)(» -  <■) 


««-  <« 

»-• 


SIB  ROUTINE  VALUE  <  X, Y* T  ,  4L  ,BE  (U.V,  kwmmm) 

THIS  IS  tARLE  ENTRY  NUMBER  16  IF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REAL+8  < A,d,C-H,0-Z) 

IMPLICIT  COMPLEX* 16  (Cl 
MPMMM  -  MMmmm  *  I 
ZERO  >  0.0*0 


CS1  «  OCMPLX(X.Y) 

C$2  »  DCMPLXJ BE, ZERO) 
C  S3  ■  OCMPLXIAL.ZERO) 


:Sexp(X^x»t  i 


XPTR, X0TI 1 


XPTR  -  X*T 
xri  *  Y*T 
CXPT  -  DC 
CMUM  »  cc --  .  ... 

CS4  -  CS1  -  CS3 
CS5  «  CSI  -  CS2 
COENOM  .  CS4*CS  5 
CRES1.T  -  CN’JM/COE NOM 
U  -  JREAL(CPESLT) 

V  *  DIMAG(CRESLT) 

RETURN 

END 


AL  ■  1.25D-1 
BE  -  -1.25D-1 
T  *  8.D+0 
AA  -  0,4 
Xy  -  -4.89 
yF  -  2.3 
MMMMM  -  24 
LOSC  *  4 

SUM  »  9.401610*0 
SUMA  «  1.07718D+1 
TIME  •  0.766480-1 
ANALYTICAL  •  9.401610+0 


I 

Case  17. 


•  +  * 


»«**  -  «c“ 

»-« 


SUBROUT  INE  V  ALJ  E < X , Y , 7 , 4L ,  3E , J  ,V  ,  mmm«M  ) 

This  is  table  entry  number  w  of  laplace  transforms  by  spiegel 

MPLICIT  REALMS  < A,3,0-H,0-Z) 


I MPLI Cl T  C0MPL£X*X6  (CJ 
MMMMM  ■  MMMMM  ♦  l 
ZERO  »  O.O+O 
CSl  -  OCMPLX(X.Y) 

CS2  «  OCMPLXl  BE  .ZERO) 

C$3  ■  OCMPLX (AL, ZERT) 
XPTR  »  X*T 
XPTI  -  Y*T 

CXPT  ■  OCMPLX (X°T  R»  XPTI ) 
CNUM  ■  CSl* COEXP ( CXPT) 

Els4 :  Eli :  VA 

CDENOM  «  CS4*CS5 
CRESLT  *  C'JUM/COENOM 
U  ■  OR E AL( CPESLT) 

V  «  OIMAG(CRESLT) 

RETURN 

END 


AL  -  1.25D-1 
BE  *  -1.25D-1 
T  ■  8.D+0 
AA  ■  0.3 
Xy  ■  -*.99 
yF  •  2.3 

MMMMM  -  2* 

LOSC  *  5 
SUM  ■  i.5*308D+0 
SUMA  -  2.52U01D+0 
TIME  •  0.78988D-1 
ANALYTICAL  •  1.5*308D+0 
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Case  18 


»fn«<  —  at  cm<t< 
2t> 


SUBR  OUTI NE  VALUE  (X ,Y ,T , ALtU, V , MMMMM) 

THIS  IS  TABLE  ENTRY  NUMBER  18  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REAL *8  (  A, 3 ,  C-h, Q-Z  ) 

IMPLICIT  C0MPlEX*16  (Cl 

4MMMM 


IMPLICIT  COMPLEX+16 

MMMMM  -  MMMMM  +  1 

ZERO  -  O.O+O 
CSI  ■  0CMPLX(X,Y1 
XPTR  *  X*T 
XPTI  •  Y*T 

CXPT  •  DCMPLX(XPTR.XPTI) 
CNUM  *  CDEXP(CXP7  1 
C$4  -  CS1*CSI 
CS5  •  DCMPLXIAL.ZEROl 
CS6  ■  CS5+CS5 
C$7  -  CS4  ♦  CS6 
COENO*  -  CS7MCS7 
CRESLT  «  CNUM/CDEMOM 
U  »  OREAL(CRESLT) 

V  »  0IMAGICRESLT1 
RETURN 


AL  «  1.25D-1 
T  »  8.0+0 
AA  -  0.3 
Xr  ■  -4.11 
yp  »  2.1 
MMMMM  -  22 
LOSC  *  2 

SUM  ■  7.70992D+1 
SUMA  «  7.70992D1 
TIME  •  0.51870D-1 
ANALYTICAL  7 . 709920+1 
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?HlS°i'sIT4BLEL6NTRY,N6UER,19V6rrAPL4CE  TRANSFORMS  BY  SPIEGEL 


IMPLICIT  REAL  *8  (  A,  B,  O-ti,  O-l  I 
IMPLICIT  COMPLEX*  16  <  Cl 
NMMMM  «  MMMMM  ♦  1 

ZERO  »  0.0*0 
CS1  -  OCMPLX(X,rj 
XPTR  »  X*T 
XPTI  *  Y* 

CXPT  ■  OCMPLXIXPTR.XPTI) 

CS3  «  COEX°(CXPT) 

;njm  •  csi*c S3 

54  *  C  S1*C  SI 

55  *  DCMPLX(AL.ZERO) 

56  *  CS5*CS5 
C$7  «  CS*  +  CS6 

SSigPT : 

V  :  DiSAbStRiltr I 

RETURN 
END 


AL  m  1.25D-1 
T  ■  8.D*0 
AA  ■  0.3 

XF  -  -4.54 
yy  »  2.2 
MMMMM  -  23 
LOSC  *  3 

SUM  *  2.692710*1 
SUMA  •  2.69278D+1 
TIME  ■  0.89190D-1 
ANALYTICAL  2.69271D*! 


Case  20 


■1 


»> 


lin  •<  +  at  cot  et 

2m 


?H1  S°Vs ^  T  f  8 LE  *"E Nt R  Y *  NU M9E R '  2 6^ 0  F **  L  A®  LAC E  TRANSFORMS  BY  SPIEGEL 


mmmmm 

ZERO 


M'4M  MM 
0.0*0 


SI  •  OCMPLXIX.Y) 

CS2  ■  CS1*CS1 
XPTR  »  X*t 

CXP^  •  DCMPLXIXPTR.XPTI ) 
tS3  «  COEX P(CXPT ) 

JNUN  -  CS2*CS3 
”>4  ■  CS1*TS1 
iS  *  0CMPLX(AL«2ER0) 

6  •  CS5*CS5 
_J7  »  CS4  ♦  CS6 
COE  NON  -  CS7*C$7 
CRESLT  «  CNUM/COENOM 
U  «  DREAL  (CRESLT) 

V  *  01 MAG(CRESLT) 

RETURN 

END 


AL  -  1.25D-1 
T  »  8.D+0 
AA  -  0.3 
XF  -  -4.54 

yp  •  2.2 

MMMMM  •  23 
LOSC  »  4 

SUM  ■  5.52709D+0 
SUMA  ■  5.89548D+0 
TIME  ■  0.60684D-1 

analytical  5.52709D+0 
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Case  21 


-1 


<«*+•*>« 


eo ««/  —  *in  at 


SUBROUT IME  VALUE  ( X, Y, T. 4L , U, V, MMMMK) 
THIS  IS  TABLE  ENT RY  NUMBER  A  OF  LAPUCE 
IMPLICIT  REAL*8  ( 4,3 .O-H.O-Z) 

IMPLICIT  COMPL  EX* 16  (C ) 

MMMMM  »  M«MMM  ♦  I 
ZERO  «  0.0*0 
CSi  -  DCMPLX(X,Y) 

CS1  -  OCMPLX(X.Y) 

CS2  *  CSl*CSl*C5i 
XPTR  -  X*T 
X  PTI  *  Y*T 

CXPT  -  DCMPLXI  XPTR,  XPTI  ) 

CS3  »  CDEXP(CXPT) 

CNUM  -  CS2*CS3 

CS4  »  CS1*CSI 

CS5  «  DCMPLXIAL, ZERO) 

C  S6  *  C  S5*C  S5 
CS7  «  CS4  ♦  CS6 

C?fsLT  «  ^NUM/!ogNOM 
U  ■  DREAL(CRESLT) 

V  -  DIMAGICRESLT) 

RETURN 

END 


TRANSFORMS  BY  SPIEGEL 


AL  -  1 .250- 1 
T  -  8.D+0 
AA  -  0.3 
Xp  -  -4.99 
y7  -  2.3 
MMMMM  -  24 
LOSC  »  5 

SUM  ■  1.19567D-1 
SUMA  •  1.S9481D+0 
TIME  •  0.10865D+0 
ANALYTICAL  1.19567D-1 
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Case  22 


SU8R"MJT INE  VALUE  IXtYfT.ALtU.V.MMMMM) 

THIS  IS  TABLE  ENTRY  NUMBER  22  jF  LAPLACE  TRANSFORMS  BY 
impli:it  RE 4L*8  ( A.B.O-H.O-Z) 

I “PL I C IT  C1M0LEX*16  1C) 


SPIEGEL 


MMMMM  ♦  1 

0.0*0 


MKMMM 

ZERO  ■  - 

CS1  •  OCMPLXIX.Y) 

C«2  -  OCMP LX( AL t ZERO) 

XPTR  ■  X*T 
XPTI  «  Y*T 

CXPT  *  OCMOLXIXPTR.XPTI) 

CNUM  -  cTcsWsi)  -  (CS2*CS2>  J*CO£XP(CXPT) 
CS3  -  MCSl*CSl)  ♦  (  CS2*CS2 ) ) 

CDENOM  -  C$3*CS3 
CRESLT  «  CNUM/C06N0M 
U  ■  ORE AL(CR£SLT ) 

V  •  01 M AG( CRESLT) 

RETURN 

ENO 


AL  -  1.25D-1 

T  -  8.0+0 

AA  »  0.3 

Xp  -  -4.51+ 

yp  -  2.2 

MMMMM  -  23 

LOSC  »  4 

SUM  •  4. 32242D+0 

SURA  -  4.76000D+0 

TIME  »  0.67210D-1 

ANALYTICAL  4.32242D+0 
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Case  23 


«t  eothtt  —  tinhof 

1? 


IMPLICIT  *  E4L *8  ( A,3,0-H,0-Z) 
IMPLICIT  COMPLEX "*16  (C) 

MMMMH  -  MMMMM  *  l 
ZEPO  •  0.0*0 
Cfl  •  ocmplxix.y) 

XPT»  >  XPT 
ert  m  y  *7 

DCMPLXIXOTR.XPTI J 
;tJE  XP(CXPT) 


-  C$l*CSi 
•  OCMPLXIAL.ZcPO) 
;S6  «  CS5*CS5 
J$7  -  CS4  -  CS6 
CDENOM  -  C  S7*C  S7 
CRESLT  -  CNUM/C0EM3M 
U  »  0 RE AL (C PESLT) 

V  -  OIMASICRESLT) 

RETJRN 

END 


AL  »  1.2 50-1 
T  ■  8.D*0 
AA  »  0.4 
XP  -  -4.44 

yF  -  2.2 

MMMMM  »  23 

LOSC  «  2 

SUM  «  9.417710*1 
SUMA  ■  9«i»1771D*1 
TIME  ■  0.61984D-1 
ANALYTICAL  9.41771D*! 
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Case  24 


<  linh  af 


U 


teansfo.«  er  Siegel 

IKPUCIT  REAi.*§  (A,3.C-H,>Z) 


jMSLICIT  C0MP|_EX*16  (C) 
MPMMM  «  ♦  1 

ZERO  -  O.O+O 
CSl  •  DCMPLXIX.Y) 

X  PTR  »  X*T 
YD  T!  m  Y*t 

';molx<xptr.xpti ) 


XOtCXP"') 


Eirc-sc5s 


ZERO) 


CS6 

JOENO't  •  CS  7*CS7 
CRESLT  -  CNUM/COcNOH 
U  «  oreal<creslt; 

V  -  01  MAG JCRESLT) 
R6TJRN 
6  NO 


AL  »  1.250-1 
T  ■  8.D+0 
AA  ■  0.4 
Xp  -  .4.U4 

yp  ■  2.2 

MMMMM  ■  23 

LOSC  ■  3 

SUH  -  3.76064D+1 
SUMA  •  3.76151D+1 
TIME  •  0.66924D.1 
ANALYTICAL  3.76064D+1 
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Case  25 


tilth  at  +  at  cosh  at 
U 


»y  sp.£ki 

1MPLI7IT  R£AL*8  ( A, B, 0-H, 3-Z > 


1PPLKIT  COMPLEX*  16  (C) 
MPMMM  ■  M*MMM  ♦  1 
ZERO  *  0.0*0 

csi  ■  0CMPtx(x,y) 

CS2  *  CSl*CSl 
XPTR^«  X*~ 

CXp4  -  OCMPLXCXPTR.XPTI) 
CS3  *  COEX°( C XP T ) 

CNUM  ■  CS2*CS3 
CS4  »  CS1*CS1 
CS5  »  OCM°  LX  <  AL  *Z  ERO ) 

ccI6t  : 

CRESL*  *  ^NUM/COENOM 
U  «  3REALCCRESLT) 

V  *  OIMAGICRESLT) 

RETURN 
END 


AL  -  1.25D-1 
T  -  8.D*0 
AA  •  0.4 
Xp  -  -4.89 

yp  -  2.3 

MMMMM  ■  24 

LOSC  -  4 

SUM  ■  1.08731D+1 

SUMA  •  1.208700*1 

TIME  -  0.73710D-1 

ANALYTICAL  *  1.O0731IH1 
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-  '  .  M...*,  , 


Case  26 


cosh  •«  +  |a<  (inh  at 


imhsfoiws  >.  »i«el 

IMPLICIT  REAL*8  U,B.C-H,W) 

illicit  complex* i6  (ci 

mhmmm  »  MMMMM  ♦  1 
ZERO  «  0.0*0 
CS1  »  OCMPLXI  X.  Y) 

CS2  *  CSI*CS1*C$1 
XPTR  *  X*T 
XPTI  *  Y*T 

CXPT  ■  OCMPLXIXRT  R.XPTI) 

C S3  «  CDEXP(CXPT) 

CNUM  »  CS2*CS3 
CS4  «  CSI*CS1 
C$5  *  DCMpLX( AL  »  ZERO  J 
C$6  *  C$5*CS5 
CS7  «  C$4  -  C$6 
COEM3M  ■  C  57*C  S  7 
CRESLT  -  CNUM/COEM1M 
U  -  ORE ALCCRESLT) 

V  ■  0 1  MAGIC  RE  $LT) 

RETURN 

EM) 


AL  -  1.25D-1 
T  *  8.D+0 
AA  «  0.4 
Xp  *  -4.89 
yr  -  2.3 

MMMMM  •  24 

Lose  *  5 

SUM  ■  2.13068D+0 
SUMA  •  4.30908D+0 
TIME  ■  0.19^90D+0 
ANALYTICAL  2.13068D+0 
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Case  27 


$t  *  a- 


tcotkat 


SUBROUTINE  VALUE  (  X,  Y,  T, AL, U,  V .MMmvmj 

THIS  IS  TABLE  ETTRY  NUMBER  27  OF  LA°LACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  R  6  AL  *  8  { A, 8 , 0-F ,0-Z ) 

IMPLICIT  COMPLEX* 16  (CJ 
MMMMM  «  MMNMM  ♦  1 
ZERO  «  0.0*0 
CS1  *  DCM»LX(X,  Y) 

CS2  *  OCMPLXIAL.ZERT) 

XPTR  *  X*T 

2xpf  -  ?CMDLXIXPTR,XPTI) 

CNIM  ■  <CCSl*CSl)  +  (CS2*CS2))*CDEXP(CXPT) 

CS3  *  UCSl*C$i)  -  <  CS2*CS2) ) 

COENOM  ■  C  S3*CS  3 
CRESLT  »  CNUM/COcNOM 
U  »  OR EAL ( rR ESLT ) 

V  -  OIMAG(CRESLT) 

RETURN 

END 


AL  «  1.25D-1 
T  ■  8.D+0 
AA  -  0.4 
XF  -  -4.89 
yF  *  2*3 

MMMMM  -  24 

LOSC  ■  4 

SUM  «  1. 23446 D*1 

SUMA  «  1.34023D+1 

TIME  -  0.65078D-1 

ANALYTICAL  1. 234460*1 


19 


THIS  IS  ’ABLE  ENTRY  NUMSI 
I MPLIC IT  PEAL *3  <A,8»C-h 
IMPLICIT  COMPLEX* 16  (C) 
HPMMM  m  MMMMM  ♦  l 

ZERO  *  0.0+0 

csi  *  ocmplxix.y) 

XPTR  «  X*I 
v  PTT  *  Y*r 

CXPT  «  DCMPLXIXPTR.XPTI) 
CKJM  «  COEXP(CXPT) 

c Is  *  OCMPLXI AL  tZERO ) 

C  S7  :  CCir+CS^S6 

COE*n«  *  CS7*CS7*CS7 
CRESLT  «  CNUM/COENOM 
U  -ORE  AH  CRESLT) 

V  -  DIMAGICRESLT) 

RETURN 

END 


LOSC  »•  2 

SUM  *  2.5^0960+2 
SUMA  *  2.89148D+2 
TIME  *  0.50U14D-1 
ANALYTICAL  «  2.54096D+2 


Case  29. 


(«-  +  «’)■’ 


f  tin  at  —  ttt-  cos  at 
8r» 


SUBROUTINE  VALUE  (X  , Y  ,T  ,  AL,U ,  V  ,  MMMf/M ) 

THIS  IS  TA  OLE  ENTRY  NUMBER  29  OF  LAP LACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REAL *3  < A,8,D-F,7-Z) 

IMPLICIT  COMPLEX*  16  (Cl 

MMMMM  <  MSMMM  ■*.  1 

ZERO  «  O.O+O 
CSI  «  OCMPLX(X,YJ 
XPTR  ■  X*T 
XPTI  *  Y*t 

CXPT  «  DCMPLXCXPTR.XPTI J 
CS3  ■  COEXO(CXPT) 

CNUM  ■  CS1+CS3 

CS4  *  CS1+CS1 

CSS  «  OCMPLXIAL, ZERO) 

CS6  -  CS5+CS5 
CS7  *  CS4  ♦  CS6 
CDEN3M  -  CS  7*CS7*CS7 
CRESLT  ■  CNUM /COE  NOM 
U  ■  0°  EAL ( CP  ESL  T ) 

V  *  DIM AG( CRESL T I 

RETURN 

END 


AL  -  1.25D-1 
T  -  8.D+0 
AA  ■  0.3 
Xp  -  -4.11 

yF  -  2.1 

MMMMM  ■  22 
LOSC  *  3 

SUM  »  1.54198D+2 
SUMA  ■  1.66758D+2 
TIME  *  0.56472D-1 
ANALYTICAL  •  1.54198D+2 


Case  30. 


<£'x 


(t2  +  ••)» 


” *  (1  +  tin  a.  -  at  cat  a! 

-  8#1 


‘-Isl  el^t  i  ^  v6FMr:si4CE  i»«SFm« «. 

IMPLICIT  REAl*8  l A , B , O-H ,0-Z) 

IMPLICI"  - - - 


JMPLICIT  C0mPL£X*16  <C) 
MMMMM  *  MMMMM  «■  1 
ZERO  -  O.CH-O 
csi  «  oc**PLXtx,r> 

CS2  «  CS1*CS1 
XPTR  «  X*T 
XOTI  *  Y*T 

CXPT  ■  DCMPLX  (XPTRtXPY  I ) 
C  S3  -  CDEXP(CXPT) 

CNUM  «  CS2*CS3 

CS4  -  CS1*CS1 

C|5  •  DC^PLXUL  .ZERO) 

CS7  «  CS4  ♦  CS6 
COENPM  »  CS7*CS7*CS7 
CRESLT  -  CNUM/COENOM 
U  *  OREALICRESLT) 

V  »  DIM AG(CRESLT) 

RETURN 

END 


AL  «  1.25D-1 
T  -  8.D+0 
AA  •  0.3 
XF  *  -^.11 

yF  -  2-1 

MMMMM  -  22 
LOSC  -  2 

SOM  ■  7.312890*1 
SUMA  -  7.31289D+1 
TIME  -  0.57122D-1 
ANALYTICAL  -  7.31289D+1 
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Case  31 


31  sin  mt  1-  at:  ?o*ut 

”  Urn 


SUBROUTINE  VALUE  t X, Y, T, At, U, 

This  JS  TABLE  E'lTRY  NUMBER  31  OF  LAPLACE 
IMPLICIT  5?(U8  t  A,9 ,3-H.O-Z) 

IMPLICIT  C0Mt>LEX*l6  (Cl 
MMMMM  «  MMMMM  ♦  I 
ZERO  «  0.D+0 
CS1  »  0CMOLX ( X.  Y) 

C$2  »  C$1*CSI*CS1 
X°TR  *  X*T 
XPTI  «  V*” 

CXPT  •  0CMPLX(XPTR,XPT1I 
C S3  »  COEXR( C  X°  T ) 

CNUM  ■  CS2  *CS3 

C$4  -  C$l*CSl 

C$5  *  OC^LXUUZERO) 

C$6  ■  CS5*C$5 
C$7  -  CS4  ♦  CS6 
C DEMON  »  C*  7*CS  7*CS7 
CUESL"  *  CNUM/Ci>fcNOM 
U  «  3«EAL(CRESLt) 

V  -  OIMAG(CRESLT) 

RETURN 

END 


TRANSFORMS  BY  SPIEGEL 


AL  »  1.25D-1 
T  •  8.D+0 
AA  »  0.3 
XF  -  -4.54 
yF  -  2.2 
IHMKMM  *  23 
LOSC  «  3 

SUM  «  2.U5177DT1 
SUMA  •  2.45184D+1 
TIMS  »  0.60554D-1 
ANALYTICAL  «  2.451770+1 
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Case  32 


7MMF0MS  .»  SPI EGEL 


I MPLIC lT  fl £ AL*8  < A.3.C-h,0-Z) 
IMPLICIT  COMPLEX’*  16  (C) 

H*MMM  ■  HMMMM  >  1 

ZERO  ■  0.0*0 
CS1  ■  OCMPLXU.Yl 

ifinVUfc*1*6* 1  *cs  1 

XPTI  »  V*T 

CXPT  ■  OC^PLX  (XPT  R»XPT  (  ) 

C$5  ■  OcipLX(  AL.ZEROJ 
CS6  ■  C$5*TS5 
CS7  •  CS4  ♦  ■  CS6 
*06N0M  -  CS7*CS7*CS7 
CRESL T  *  CVUM/COSNON 
U  -  DREAL(CRESLX> 

V  -  DINAGICRESLTJ 

RETJRN 

END 


AL  »  1.25D-1 
T  •  8.D+0 
AA  •  0.3 
Xp  •  -4.54 

yp  •  2.2 

MMMMM  «  23 

LOSC  -  4 

SUM  •  4.38445D+0 
SUMA  ■  4.83048D+0 
TIME  ■  6.7548D-2 
ANALYTICAL  »  4.38445IH0 


SUBROUTINE  VALUE  (X »Y »T , ALt J *V . mmvmm) 
THIS  IS  TA3LE  ENTRY  NUMBER  33  OF  LAPLACE 
IMPLICIT  ft  E AL*8  ( A, 9i O-H.O-Z ) 

I  PPUC I  T  .COMPLEX*  16  (  C ) 


TRANSFORMS  BY  SPIEGEL 


MMM 


m  mmmmm  ♦  i 


ZERO  ■  0.0*0 

ey :  ?ir*hsxi*xd?i*c 

X  PTR  -  X  *T 


l*CSl*CSl 


NON 


.  Y*T 

»  OCMPLXl  XPTR,  XPTI ) 

COEX  0 (CXP“ ) 

•  C  S2*C  S3 

OCMPLX ( AL. ZERO! 

9rs,~ 


Islt  -  SNUM/l 

•  SimagIErII 


RETJPN 
END 


AL  -  1.25D-1 
T  *  8.D+0 
AA  *  0.3 
Xp  -  -4.99 
yF  -  2.3 

MMMMM  *  24 
LOSC  »  5 

SUM  -  -2.63523D-1 
SUMA  ■  1.81943D+0 
TIME  •  6.74414D-2 
ANALYTICAL  -  -2.63523D-1 


Case  3^ 


3a*  —  a* 

<j* +«j? 


t*  »ln«< 

u 


?«?|0ySI?fB^L6§TW,Nili5ML,3;''&rn?iACE  TRANSFORMS  BV  SPIEGEL 
IMPLICIT  REAL *8  ( A, 9. 0-H, 0-Z ) 

IMPLICIT  C0MPLEX*16  CC> 

mMMMM  w  MMMMM  ♦  1 

ZERO  *  0.0*0 

m  :  »pllx*IaxcUero» 

XPTR  »  X*Z 
XPTI  «  Y*T 

CXPT  -  pCMPLXIXPTR,XPTII, 

Ss3*-  . 

CrIslT  »  CNijM/CDENOM 

U  ■  0°  EAL(CRESLT) 

V  -  OIMAG(dRESLT) 

R  ETUR  N 
ENO 


R  •  A  ^ 

I  *  Y*T 

\  i  « 5:W45?  11  li?|II,-,cs2.rsi»»-coex.ccxPT. 
-  !CS1*CS1)  *  (CS2»CS2) 

NOM  ■  CS  3*CS  3*CS  3 


A!  -  1.25D-1 
T  «  8.D+0 
AA  -  0.4 
Xy  «  -4.44 

yp  ■  2*2 

MMMMM  ■  23 

LOSC  -  2 

SUM  »  2.154170*2 
SUMA  ■  2.154170*2 
TIME  -  0.678340-1 
ANALYTICAL  *  2.15147D+2 
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Case  35. 


«»-  3aU  - 

(PTZ?  — 


j/1  coaat 


TRANSFORMS  .Y  SP  I EGEL 

IMPLICI  T  R£4L*8  <  4,0  ,0-H,0-Z) 

IMPLICIT  COMPLEX* 16  (C) 

MMMMM  *  MMMMM  ♦  1 
ZERO  »  0.0+0 
CSI  «  DCMOLXtX.r) 

CS2  •  DCMPLXIALt ZERO) 

XPR  -  X*T 

£nUM  •  ??C S 1*CS l*c£ lf-( 3. D+0*CS  t*CS2*C  S2 ) ) *C0EXP (CXPT ) 

CS3  -  <CS1*CS1)  ♦  ( CS2*CS2) 

COENOM  .  CS  3*CSi*  CS3 
CRESLT  »  CMUM/C3EN0M 
'J  -  ORE  AL  ( CPE$LT ) 

V  *  DIMAG(CRESLT) 

RETURM 

END 


AL  -  1.25D-1 
T  ■  8.D+0 
AA  *  0.3 
Xr  •  -4.54 
yp  •  2.2 
MMMMM  *  23 
LOSC  -  3 

SUM  ■  1.72897D+1 
SUMA  *  1.72904D+1 
TIME  •  8.44740-2 
ANALYTICAL  ■  1.72897D*! 


Case  36. 


t*  -  w  *  «« 


£(•«»«< 


hmmnm ) 

f  LAPLACE  TRANSFORMS  BY  SPIEGEL 


MPMMM  •  MMUNM  *  I 
ZERO  »  0.0+0 
CS1  -  DCMPLXIX.Y) 
CS2  ■  OCMPLXIAL.I 
C S3  -  CSI*CS1*CS1 
C$4  -  6.D+OK5Z*: 
CSS  »  CS2*CS2*CS. 
XPTR  »  X*T 


OC  MPLXI X»  Y) 
OCMPLX (AL.ZE90 

csi*csi*cn*c« 

6. D+O+C  52*C  S2* 


2*CS2*CS2*CS2 


Sl*CSl 


Cxpf  »  OCMPLX (XPTR.XPT I J 

CNIM  «  ICS3-CS4+C S5)*CDEXP(CXPTJ 
CS6  -  <(0S1*CS1>  +  ( C  52+CS2)  ) 
CDENOM  ■  CS6*C5o*  CTo*CS6 
CRESLT  ■  C  NUM/COE  NON 
U  ■  ORE AL (C  RESLT ) 

V  »  DI MAGICRESLT) 

RETl*N 
END 


OENOM  ■  CS6*C5o*  CTo*CS6 
*ESLT  ■  C  NUM/COE  NON 
■  OREAL (CRESLT ) 


AL  •  1.25D-1 
T  ■  8.0+0 
AA  »0.4 


yp  *  2.2 

MMMMM  *  23 

LOSC  *  2 

SUM  -  4.610580+1 
SUMA  *  4.610580+1 
TIME  -  0.73424D-1 
ANALYTICAL  -  4.61058D+1 


■  »  JP*.'.-  *»PWWWW 


Case  37 • 


H  tin  at 
24« 


SUBROUTINE  VALUE  { X .  Y  , T  ,AL.U,  V  .MMMMM) 

▼HIS  IS  TABLE  E  47  AY  NU 43ER  37  OF  LAPLACE  TRANSFORMS  BY  S®I  EGEL 
IMPLICIT  R E AL*8  l A,3 ,0-H,3-2) 

IMPLICIT  00  4°LEX* 16  (C) 

MMMMM  •  MMMMM  l 
ZERO  «  C.0*0 
CS1  «  DCMPLX(X.Y) 

C^«a0CMPLX(AL.ZER0) 

XPTI  m  y*T 

CXPT  -  0CMPLX(X»T«,XP’-1) 

i>SV(Ui!i;Sit'55Hsi>Si!f*c”*tsl"'c“E  ,a  ’ 

COENOM  -  C  S3*CS3*  CS3*CS3 
CRESLT  -  CNJM/COENOM 

uv :  8?m  xbJ  ^rIIltI 

RETURN 

END 


AL  ■  1.25D-1 

T  ■  8.IH0 

AA  ■  0.4 

Xy  ■  -4.01 

yF  -  2.1 

MMMMM  ■  22 

Lose  ■  3 

SUM  •  1.43611D+2 

SUMA  *1.44005D+2 

TIME  ■  0.68692D-1 

ANALYTICAL  -  1.43611D+2 


Case  38 


(3  +  «»t*)  linh  at  -  3«<  folh  at 

s? 


BY  SPIEGEL 

i'lPLICT - -  *  *  “  ~  ' 


.•ILLICIT  f  E4L*8  <  A, 8 
IMPLICIT  C0‘13L EX* lo 
M'WMM  *  MMM*M  •*•  1 
ZEPT  *  0.0*0 
CS1  *  0CMOLXCX.Y1 
xpn  »  x*t 
XOTI  «  Y*T 

CXPT  ■  CC-tPLX  (XPTP.XOT  I) 

cnom  «  coex»<cxpn 

CSA  »  CS1*CS  1 

C$5  «  DCMPUUAL.ZERO) 

CS6  «  CS5*CS5 
CS7  »  CS4  -  CS6 
COENO*  -  CS 7*CS7*CS7 
C«6SLT  ■  rviuM/COEi‘10'1 
U  •  OPEAL(CP.ESLT) 

V  «  DI  HAG(CRE3Lt) 

RETJ9N 

ENO 


t  0_H»0-Z ) 
I  Cl 


AL  -  1.25P-1 
T  ■  8.D+0 
AA  ■  0.4 
Xp  -  >4.01 
yp  -  2-i 

MMMMM  »  22 
LOSC  -  2 

SUM  -  2.93122D+2 
SUMA  -  3.67471D+2 
TIME  -  0.48334D-1 
ANALYTICAL  *  2.93122D+2 
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Case  39 


«l»co«h«(  ~  t  linhat 

la* 


SUBROUTINE  VALUE  (X,  Y.T.ALtU,  V.VMMMM) 

THIS  IS  ’’ABLE  ENT  ST  NUMBER  39  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REALMS  (  A, 9 . Q-H ,0-Z) 

I  PPL  I C IT  CnN0LEx*16  tC> 

MPMMM  »  MNMMM  ♦  l 
ZERO  «  C.D+O 
CSl  -  OCM°LX (X*  Y I 
XPTR  «  X*T 
XPTI  »  Y*T 

CXPT  *  OCM°LXCXPTRfXt'Y|J 

CS3  -  COEXP( C XPTI 

CNUM  ■  CSI+CS3 

CS4  •  CS1*CS1 

CSS  ■  DCMPLXIALtZEROl 

CS6  »  CS5*CR  5 

C$7  ■  C$4  -  CS6 

COENOM  -  C$7*CS7*C$7 

CRESL’-  «  CNUM/CJENOM 

U  »  OREAL(CRESlT) 

V  -  OIMAGICRESLTJ 

RET,JRN 

END 


AL  -  1.25D-1 
T  •  S.D+O 
AA  »  0.4 
Xp  -  .4.01 

yP  •  2.1 

MMMMM  ■  22 
LOSC  -  3 

SUM  ■  1.88354D+2 
SUMA  •  1.89195D+2 
TIME  ■  0.48750D-1 
ANALYTICAL  ■  1.88354D+2 


131 


Case  40 


-1 


«« 


mt  eu*b  *t  +  (a*/*—  1)  linh  it 
8a> 


c 


ftfiriSI!r!pLEL^TS^M^d^^fMU?CACE  TRANSFORMS  BY  SPIESEL 
-  -CIT  REAL*B  ( A.BfC^H.O-Z) 


IMPLI  n 
IMPLI-IT  COMPLEX*  16  4  c> 
MRMMM  »  MMM“M  ♦  1 

ZERO  ■  O.O+O 
CSl  *  OCMPLXI  X»  Y) 

CS2  *  CS1*CS1 

xptr  -  x*i 

XPT1  *  Y*T 

CXPT  -  0CM3LX(XPTR,XPTI) 

csa-  coix;.c«T, 


CS1*CS1 
0CMPLX(AL.ZER01 


CS4 

ib :  aeas 


C  0EN3N 
CPESL 
u  «  o 

V  -  OIMA 
RETURN 
END 


C  *7*C  $7*C  $7 
r.NUM/COENOM 
,  (CRE5LT ) 
i  ( C  a  E SLT) 


AL  -  1.25D-1 
T  »  8.D+0 
AA  •  0.4 
X?  -  -4.44 

yp  -  2.2 
MMMMM  *  23 
LOSC  -  2 

SUM  «  9.875720+1 
SUMA  *  9.875720+1 
TIME  *  0.52520-1 
ANALYTICAL  *  9.87572D+1 
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T..NSFW-S  .Y  SPIEGEL 


mmmmm  ■  MVMMM  +  i 

ZERO  =  0.0+3 

CS1  *  DCMPLX(  X.Y)  , 

CS2  *  CSl*CSl*CSi*CSl 
XPTR  »  X*- 

CXPT  *  0c5°tx  (XOTR.XPTH 
CS3  »  COEXP(CXPT) 

GNU*  -  CS2*CS3 

Els  «  0CMPLX( AL.ZERO) 

CS6  =  CSS+'-SS 
CS7  -  CS4  -  C36 
CDENDM  -  C  S7*C  S7*CS7 
CRESLT  »  CNUM/CDr  0* 
u  »  ORE  AL ( CRESLT  t  . 

V  »  OIMAG(CRESLT) 

RETURN 

END 


AL  -  1.25D-1 
T  »  8.D+0 
AA  -  0.4 
Xp  -  -4.89 
Yf  «  2.3 

MMMMM  »  24 
LOSC  -  4 

SUM  »  1.24162D+1 


SUMA  -  1.34824D+1 


T.A..SFOB..S  BY  SPIEGEL 

implicit  peal *8  < a* b. c-Eto-z ) 
implicit  cn«PLex*ifc  to 

nmMMM  a  1 

ZERO  ■  0.0+0 

c§2  *  csi*csi*dsi*csi*csi 

XPTR  «  X*T 

ffl  »  CCMPLXIXPTR  ,XPTI  1 
CS3  »  CDEX °( CXPT ) 

CNUR  »  C S2+CS3 

£§5  «  DCMPOX ( ALtZERO) 

C$6  *  CS5*CS5 
CS7  ■  CS4  -  rSb 

BKSe? ;  5®55?sa? 

if  -  ORE  At  ICRES17 ) 

V  »  OIMAC(CRESLT) 

RETURN 

END 


AL  *  1.2 50-1 
T  »  8.D+0 
AA  *  0.4 
Xjp  ■  -4.89 
yF  -  2.3 

MMKMM  -  24 

LOSC  ■  5 

SUM  -  2.7^4270+0 
SUMA  -  4.73261D+0 
TIME  «  0.629200-1 
ANALYTICAL  -  2.?6427D+0 
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Case  44. 


>«w»i  •» 

IMPLICIT  e  E^L *3  ( A, 6 ,0-H,0-Z ) 

IMPLICIT  C0*<PLcX*l6  ( Ci 

UM*MM  a  *  l 

ZE«T  *  0.0+0 
C51  «  0CMPUX(XfY) 

CS2  «  0CMPLX(AL,Z6R0» 

XPTR  «  X*T 
X  PTI  •  Y*T 

CUM  «  1 13^0+0* ( cl  1*CS \n+(CS2*CS2»)*C0EXP(CXPTJ 
CS3  *  CCSl+CSll  -  (CS2*CS2) 

COE'IOM  a  CS3*"S3»CS 3 
CRESLT  «  CN'Jw/COEMO* 

U  *  ORE  AL( CRE  SLT) 

V  -  OI*A&(CRESLT» 

RETURN 

END 


AL  -  1.25D-1 
T  -  8.0+0 
AA  -  0.4 
Xp  -  -4.44 

yF  -  2.2 
MMMMM  ■  23 
LOSC  “  2 

SUM  -  3«  00852D+2 
SUMA  ■  3.00852D+2 
TIME  -  0.550420-1 
ANALYTICAL  •  3.00852D+2 
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Case  45 


|l<  CMh  at 


SUBROUTINE  VALUc  <XtY*T*AL»U» 
THIS  IS  TABLE  SNTRY  N'JMBeR  45 
IMPLICIT  B  EAL*3  ( A, B. o-hfO-Z ) 
IMPLICIT  COMPLEX* lo  (Cl 


\i  mmmmm) 

OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 


lEROJ 


Implicit  complex* io 

MMM4M  *  MM^MM  ♦  1 

ZERO  *  0.0*0 
CS1  «  0CMOL X(  X»  Y) 

CS2  *  ocmplxJal.zi 

X  PTR  »  X  *T 
X  PTI  *  Y*T 

8mm  «  HCSl*CSlicinI(3.0*0*CSl*CS2*CS2n*C0EXP(CXPT) 

CS3  ■  (  CS1*C  SI)  ”  I  C  S2*C  S2 ) 

CDENOM  *  CS3*CS3*-S3 
CRESLT  *  CMUM/COcNOM 
U  «  OREAL(CRESLT) 

V  »  OIMAGICRESLT) 

RETURN 

EMO 


AL  -  1.25D-1 
T  «  8.D+0 
AA  »  0.4 
Xp  -  -4,44 
yF  -  2.2 
MMMMM  *  23 
LOSC  -  3 

SUM  ■  4.937860+1 
SUMA  ■  4.938920+1 
TIME  -  0.57772D-1 
ANALYTICAL  -  4.937860+1 
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Case  46 


j| 


1<  +  +  «« 


Jt*  coihat 


IMPLICIT  RE«L*3  <  A,3,0-H,'W) 
IMPLICIT  COMPLEX*!*  (Cl 
MMMMM  ■  *  1 

ZEa3  «  0.0*0 
CS1  ■  OCMPLXIX.Y) 

CS2  •  DCMPLXUL.ZERO) 

CS3  -  CSI*CSl*Cii*CSl, 

C$4  *  6.0*0 *CS2*C 32 *C$l*CSl 
CS5  «  CS2*CS2*CS2*CS2 
xptr  -  x*r 
XPTI  ■  Y*T 

CXPT  »  OCMPLXl XPTR, XPTI) 

CNJM  .  (CS3*CS4*CS5  I*CCEXP(  CXPT  I 
CS6  ■  UCS1*CSI)  -  (  C  ?2*C32 1  I 
CDENOM  >  CS6*CS6*C$6*C$6 
CRESLT  ■  CNijM/COENOM 
U  «  OREAL(CRESLT) 

V  -  D IM AG( CRESLT I 

RETURN 

END 


AL  •  1.25D-* 

T  -  8.D+0 
AA  -  0.4 
X p  -  -4.44 

yF  -  2.2 

MMMMM  -  23 

LOSC  *  2 

SUM  •  1.31676D+2 
SUMA  «  1.31676D+2 
TIME  -  0.65676D-1 
ANALYTICAL  -  1.31676D+2 
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Case  47 


»»  s»iemi 

ini?-£r~  -  “  ' *  * 


I MPLI Cl T  S1Sl*8  (A,9;6-HVn-Zj 
implicit  complex* 16  <ci 

M*MMM  *  MMMMM  ♦  l 
ZERO  ■  0.0*0 
c  $1  -  DCMPLXt  X »  Y) 

CS2  -  DCMPLX  <  AL  *  ZERO ) 

XPTR  «  X*T 
XPTI  «  Y*T 

CXPT  «  CC-PLX(XPTRj,XPTn 


CNUM  -  HCSlPCS^cjfri’jCSaPCSZPCSlDPCDEXPtCXPT) 


CRESLT  »  CNUM/CDcNOM 
U  -  DREALCCRESLT) 

V  -  OIMAG(CRESLT) 
RETURN 
END 


AL  -  1.250-1 
T  ■  8.D+0 
AA  -  0.4 
Xp  «  -4.01 

yP  -  2.1 

MMMMM  «  22 
LOSC  »  3 

SUM  ■  2.005680+2 
SUMA  ■  2.01796D+2 
TIME  -  0.55822D-1 
ANALYTICAL  ■  2.005680*2 
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Case  A8 


^^ISNf«UA^T^'YNUTMiEUR'U4SV5?Mr^ACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REAL *8  (  A,  R , ,0-t-.,D-Z> 

IMPLICIT  COMPLEX*I&  (C> 

MMMMM  *  RMKMM  1 

ZERO  «  0.0+0 
CS1  «  OCMPLXIX.Y) 

XPTR  *  X*T 
XP+I  m  y*T 

CXPT  *  OCNPLXIXPTR.XPTI ) 

CNUM  «  CDEXP(CXPT) 

CS4  «  CS1*CS1*CS1 
CS5  -  OCNPLXI AL.ZERO) 

CS6  *  CS5*CS5*CS5 
CDENO*  -  C34  +  CSo 
CRESLT  ■  :nUM/C36N0M 
U  «  DR  EALI  CRESLT  J 
V  -  OIMAG(CRESLT) 

RETURN 

END 


AL  «  1.25D-1 
T  »  8.D+0 
AA  «  0.4 
Xp  «  -4.44 
yp  -  2.2 

MMMMM  ■  23 
LOSC  ■  3 

SUM  *  3.14683D+1 
SUMA  -  3.14767D+1 
TIME  ■  0.48958D-1 
ANALYTICAL  «  3.14683D+2 
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Case  kS. 


«**'«  f 

■=T 


VJtt 


+  Vi  sin 


i^2i  -  /,} 


SUBROUTINE  VALUE  (X  ,Y  ,T,  AL.U.V  .MMMMMI 
THIS  IS  TABLE  ENTRY  NUMBER  49  OF  LAPLACE 
IMPLICIT  P  E  M.  *8  ( 4,3, C-l-,Wi 
IMPLICIT  C0MPLSX*16  IC) 

MMMMM  ■  MMM'tM  «■  1 

ZERO  *  O.D+O 
C  SI  ■  OCMPLX(X.Y) 

XPTR  ■  X*T 
X  PTI  *  Y*T 

CXPT  ■  OCMPLXIXPTRtXPTII 
CS3  *  COEXO(rxPT) 

cnum  -  c!i*c«  ,, 

C  S4  -  CS1*CS1*CS1 
CS5  »  OCNPLX(AL.ZERO) 

CS6  -  CS5*CS5*C$5 
C DEMON  »  C  54  ♦  C So 
CRESL”  -  CNJM/CDE  NO** 

U  -  0REALCCRE5LII 

V  ■  OIMAGICRESLT) 

nrj  RN 

ENO 


TRANSFORMS  BY  SPIEGEL 


AL  -  1.25D-1 
T  ■  8.IHO 
AA  ■  0.3 
Xp  -  -4.54 
yp  -  2.2 
MMMMM  «  23 
LOSC  -  4 

SUM  •  7.66825D+0 
SUMA  -  7.92326D+0 
TIME  -  0.52^160-1 
ANALYTICAL  »  7.66825D+0 
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Case  SO 


T.WBWWS  BY  SP 1EGCL 

IMPLICIT  RE4l*8  (  4,8  ,0-H.O-Z) 

IMPLICIT  C^plEXPIo  (C) 

MMMMM  •  MMMMM  ♦  1 

ZERO  ■  0.0>0 
CSl  »  OCMPLX(XtY) 

C$2  »  CSl*CSl 
XPTR  »  X*T 
X  PTI  ■  Y*T 

CXPT  -  OCMPLXIXPTR.XPTI) 

CS3  «  COEXPICXPTJ 
CNJM  •  CS2*CS3 
CS4  -  c$i*cst*csi 

CS5  «  OCMOLXI AL.ZER3) 

C$6  -  CS5*CS5*C$5 
CDENOM  »  'S4  f  iS6 
CRESLT  «  CM'JM/COE  NOM 
U  »  OREALICRESLT) 

V  ■  DIMAS(CRESLT) 

RETJRN 

ENO 


AL  -  1.25D-1 
T  ■  8.D+0 
AA  *  0.4 
Xp  -  -4.89 
yF  •  2.3 

MMMMM  -  24 
LOSC  •  5 

SUM  -  8.34720D-1 
SUMA  ■  3.88818D+0 
TIME  -  0.557960-1 
ANALYTICAL  ■  8.347200-1 


Case  51 


u* 


l*4*"*  —  CM 


Vi  *t 

i 


C 


SUBAOUT  INE  VALUE  (  X,  Y.  T, AL,  U ,  V , MMMMM) 
THIS  IS  "ABLE  ENT AY  NUMBER  51  OF  LA 
IMPLICIT  R£AL*8  < A,B.D-HfO-Z) 

ILLICIT  r/"40LEx*16  IC> 


PLACE  TRANSFORMS 


MMMMM 


MMNMM  +■  1 

_.D+0 
DCMPLxix.yj 

X*T 

Y*T 

OCMPLXIXOTR, XPTI  ) 
CDEXP(CXPT) 
_S1*CS1*CS1 
OCMPLX (AL. ZERO) 
CS5* 


►  CS5*CS5 

:S4  -  cs 

J*IU“/COE 
.REALICRESLT) 
OIMAG(CRESLT) 


6 

MOM 


AL  -  1.25D-1 
T  »  8.D+0 
AA  ■  0.4 
Xy  ■  -4.44 

yP  »  2-2 

MMMMM  *  23 
LOSC  »  3 

SUM  »  3.25349D+1 
SUMA  •  3. 2 5434D+1 
TIME  ■  0.61204D-1 
ANALYTICAL  -  3.25349D+1 


BY  SPIEGEL 


Case  52 


A1  ■  1.25D-1 
T  ■  8.D+0 
AA  ■  0.4 
Xp  ■  -4.89 
yF  -  2.3 
MMMMM  -  24 
LOSC  ■  4 

SUM  »  8.33492D+0 
SUMA  •  9.862830+0 
TIME  -  0.67418D-1 
ANALYTICAL  =  8.33492D+0 
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Case  53 


ti-urt  cm 


4s) 


TRANSFORM S  BY  SMEGEL. 

IMPLICIT  R6AL*8  ( 4,3,0-H,0-Z) 

IMPLICIT  C0MBLEX*1&  I  C I 

MMMMM  a  MMMM*  ♦  l 

ZERO  •  0.0  +  0 

tsi  «  DCfPlX(X,Y) 

C$2  -  CSHCS1 
XPTR  a  X*T 
XPTI  .  Y*T 

CXPT  -  0CMOLXIXPTR,XP"  II 
C S3  «  COEXPIC XPT) 

’HUM  «  CS  2*CS 3 

54  «  C$l*CCl*CSl  „„ 

55  »  DCMPLXI AL, ZERO) 

56  «  CS5*CS5*C55 
.OENOM  «  CS4  -  CS6 
;reslt  «  CWVCOEOM 
U  -  DREAL(CRESLT) 

V  -  0IMA5(CR£SLt) 

RETJ3N 

END 


AL  «  1.25D-1 
T  »  8.D+0 
AA  -  0.3 
tf  -  -*.99 

yy  -  2.3 

Mt©WM  ■  2* 

Lose  »  5 
SUM  -  1.16806D+0 
SUMA  •  2.25373D+0 
TIME  ■  0.70200D-1 
ANALYTICAL  «  1.I6806D+0 
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s...  :  ...... 


1 


jij(»in  «f  coth  at  —  eo«  «t  >inh  •«) 


SUBROUTINE  VALUE  ( X  ,  Y  ,  T ,  AL.U,  V  ,N*«MMMI 

THIS  IS  ’'ABLE  ETRY  N!J  '8ER  54  OF  L >BL ACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REU*8  1 A,3,D~H,3-Z) 

IMPLICIT  mM(»LEX*I&  (CJ 

MHMMM  ■  MMMMM  +  l 

ZERO  «  0.0+0 

CSI  «  DCM°LX  (X»  Yi 

XPTR  »  X*T 

XPTI  -  Y*T 

CXPT  *  CCM°LX(XOTR»XOTl) 

CNUM  «  CDEXPICXPTi 


l*CSl*CSl=»CSl 

MPLXUL.ZERO) 


S4 

CS5  »  i/inr  la  |  at  {4CK  j  j 
CS6  ■  4.0+  0*( C  $5*C  $5*  CS  5*CS5) 
CDENOM  «  CS4  +  CS  6 
CR6SLT  «  CNUM/COENOM 
U  -  OREAH CRESLT) 

V  «  0 IMAGI CRESLT ) 

RETURN 

END 


AL  «  1.25D-1 
T  ■  8.D+0 
AA  -  0.4 
Xp  «  -4.44 
yF  -  2.2 

MKMMM  -  23 
LOSC  >  2 

SUM  «  8.49272D+1 
SUMA  -  8.492720+1 
TIME  *  0.677323-1 
ANALYTICAL  *  8.49272D+1 


Case  55 


SUBROUTINE  VALUE  ( X, Y, T ,AL , U, V ,MMM«M) 

THIS  IS  TABLE  ENTRY  NUMBER  55  OF  LA°LACE  TRANSFORMS  BY  SPIEGEL 
IMPLItIT  REAL*3  ( A, 3 , O-H ,0-ZJ 
IMRLICIT  timolEXmIo  (C) 

MKMMM  ■  MMMMM  ♦  L 

ZERO  *  O.O+O 
CS1  ■  OCM3LX( X» Y ) 

XPTR  »  X*T 
XPTI  •  Y*T 

CXPT  ■  DCMt>LX(XPTR,XPTI ) 

CS3  «  COEXP(CXPT) 

CNUM  -  CS1*C S3 

CS4  »  C$1*CS1*CS1*CS 1 

CSS  ■  OCMPLXI ALfZERO) 

CS6  «  4.0+0*<CS5*CS5*CS5*CS5) 

COENOM  ■  ♦  CS6 

CR6SLT  «  CNUM/CDENOM 
U  «  DREAL(CRESLT) 

V  -  DIHtGICIESLT) 

RETURN 

ENO 


AL  -  1.25D-1 
T  -  8.D+0 
AA  •  0.** 

XF  -  -A. 54 

yF  -  2.2 

MMMMM  ■  23 
LOSC  -  3 

SUM  ■  3.16447D+1 
SUMA  *  3.16532D+1 
TIME  *  0.59358D-1 
ANALYTICAL  -  3.l644?D*l 
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Case  56 


~~  (»in  at  cook  at 


+  coa  at  gink  at) 


SUBROUTINE  VALUE  ( X  ,  Y  ,T  ,  AL  ,  U,  V.MMMmm  ) 

THIS  IS  TABLE  ENTRY  NUMBER  56  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  R EAL*8  I A,B,D-H,0-Z) 

IKPLICIT  :0MPLEX*16  (C) 

MMMMM  ■  MMHMM  ♦  1 
ZERO  ■  O.O+O 
CSl  ■  OCMPLX(X.Y) 

CSZ  *  CS1+CS1 
X  PTR  «  X*’’ 

XPTI  ■  Y*T 

CXPT  -  DCMPLXI XPTR, XPTI) 

CS3  «  COEXF(CXPT) 

CNUM  -  CS2*CS3 

CS4  -  CSl*CSl*CSl*r:Sl 

CS5  -  DCMPLXI AL,ZE=OI 

CS6  *  4.D*0*ICS5*CS5*CS5*CS5) 

COENOM  «  CS4  ♦  CS6 
CRESLT  *  CNIJM/COE.NOM 
U  *  3RE  AL  ( C RE  SL  “I 
V  -  DIM AGI CRESLT J 
RETURN 
END 


AL  -  1.25D-1 
T  -  8.D+0 
AA  *  0.4 
XF  -  -4.89 
yF  -  2.3 

MMMMM  -  24 

lose  -  4 

SUM  =  7.73369D+0 

SUMA  *  9.23364D+0 

TIME  -  0.76206D-1 

ANALYTICAL  *  7.733692+0 
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Case  57 


Y3  ”i  «» 
oL,  «*  +  w 


co«  mt  coah 


SUBROUTINE  VALJE  { Xf Y, T, AL, U, V. MMMMM) 

THIS  IS  "&eLE  ENTRY  NUMBER  S7  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLI Cl T  R 5AL*8  ( 4  ,3  .0-H , 0-7) 

IMPLICIT  COMPLEX *16  (C) 

MMMMM  «  MMMMM  l 
ZERO  *  0.0+0 
CS1  *  0CMPLXIX.Y1 
C$2  »  CS1*CS1*CS1 
XPTR  «  X*T 
XPT I  «  Y*t 

CXPT  «  OCMPLX(XPTR»XPTI) 

C  S3  *  COEXP(CXPT) 

CNUM  »  CS 2*CS 3 

CS4  «  CSI*CS1*CS1*CS1 

CS5  *  DCMPLXI AL.ZERO)  . .  „  , 

CS6  *  4 .0+0* ( CS  5*  CS  5*CSa*CS5 ) 

COENOM  «  CS4  ♦  CS6 
CRESL”  ■  CNIJM/COENOM 
tl  «  ORE AL (CRESLT ) 

V  -  DIMAG(CRESLT) 

RETJRN 

END 


AL  *  1.25D-1 

T  =  8.D+0 

AA  *  0.4 

XF  -  -4.89 

yF  -  2.3 

MMMMM  ■  24 

LOSC  «  5 

SOM  *  8.33730D-1 

SUMA  -  3.85232D+0 

TIME  »  0.82576D-1 

ANALYTICAL  *  8.33730D-1 
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loot  |00003> 


Case  58 


IMPLICIT  RE  AL*8  < A,B,C-h,0-Z) 
IMPLICIT  COMPLEX* 16  (CJ 
MMMMM  m  «•  1 

ZERO  «  0.0+0 
CS1  ■  0CM°LX<X,Y> 

X  PTR  «  X*+ 

XPTI  ■  Y*T 

XPT  «  DCMPLXI XPTR.XPTI ) 

NUM  •  CDEXP( CX  ®T  ) 

54  «  CSl*CSl*CSl*C$l 

55  »  OCMPLXI  ALtZE^) 

56  «  C$5*Ci5*Ci5*CS5 
OEYOM  «  CS4  -  CS6 
RESL’  «  CNUM/COENOM 

-  OREAL(CRESLT) 

■  dimagicresltj 

RETURN 
E  NO 


AL  -  1.25D-1 
T  -  0.D+O 
AA  -  0.4 
Xy  -  -4.44 

yr  •  2.2 
MMMMM  »  23 
LOSC  -  2 

SUM  ■  8.54349D+1 
SUMA  ■  8.543490D+1 
TIME  -  0.74490D-1 
ANALYTICAL  «  8.54349D+1 


Case  59 


^j(eo«h  «<  -  cniI) 


SUBROUTINE  VALUE  (  X,  Y ,  T.AL,  U,  V  .HMMmmj 

this  is  table  emt^y  dumber  59  of  laplace  transforms  by  spiegel 

IMPLICIT  R  EAL+8  ( A, 8 . O-H ,0-2 J 
IMPLICIT  C3MOLex*16  (C) 

MMMMM  ■  MMWMM  +  i 
ZERO  »  0.0+0 
CS1  *  OCMPLXIX,  Y I 
xptr  *  x*r 

m  •  OCMPLXIXPTR.XPTIJ 
CS3  ■  COEX  P(CXPT) 

CNUM  ■  CS1+CS3 

CS4  »  CSI*CS1*CS1*CS1 

CS5  •  OCMPLX(AL.ZERO) 

CS6  «  C  S5*CS5*CS5*CS5 
CDENOM  ■  CS4  -  CSb 
CRESLT  »  CMUM/COE NON 
u  *  drealTcpesltj 

V  -  DIMAG(CRESLT) 

RETURN 

ENO 


AL  •  1.25D-1 
T  -  8.D+0 
AA  »  0.3 
Xp  -  -4.54 
*  -  2.2 
MMMMM  -  23 
LOSC  ■  3 

SUM  -  3.20889D+1 
SUMA  -  3.20897D+1 
TIME  «  0.78390D-1 
ANALYTICAL  ■  3.208890+1 
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Case  60 


1 

2a 


(linhat  +  tin  at) 


SUBROUTINE  VALUE  (X  ,Y  ,T  ,  AL  ,'J  ,  V  ) 

THIS  IS  TABLE  ENTRY  NUMBER  60  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REAL *8  ( A, A, O-F.O-Z ) 

IMPLICIT  C0MPLEX*16  (CJ 
■  NNMMM  *  i 
ZERO  *  0.0+0 
CS1  *  OCMPLXIX*  Y) 

CS2  »  CSl+CSl 
XCTR  -  X*T 
XPTI  «  Y*T 

CXPT  »  OCMPLXI XPTR t XPT I  ) 

C$3  »  COEX °ICXPT ) 

CNIM  «  CS2*CS3 

C$4  «  CSl*rSl*CSl*CSl 

CS5  -  OCMPLXIALtZERO) 

C$6  ■  CS5*CS5*CS5*CS5 
COENOM  >  CS4  -  CS6 
CRESLT  »  CNUM/COENOM 
U  -  DPEALICRESLT) 

V  ■  DIMAG(CRESLT) 

RETURN 

END 


AL  ■  1.25D-1 
T  *  8.D+0 
AA  ■  0.3 

x  -  -u.jk 
y  -  2.2 
MMMMM  -  23 
LOSC  -  4 

SUM  -  0.O6669D+O 
SUMA  ■  8.32767D+0 
TIME  »  0.87672D-1 
ANALYTICAL  -  8.06669D+0 


Case  61 


j(co«h«t  +■  oo**l) 


SUBROUTINE  VALUE  (  X,  V,  T ,  AL  ,  U,  V 

THIS  IS  TABLE  EN^RY  NU*9ER  61  OF  LA°LACE  TRANSFORMS  BY  SPIEGEL 
ILLICIT  REALMS  I  A , 8 , 0-H  ,0-2) 

IMPLICIT  COMPLEX*!*  (C) 

MPNMM  ■  MMMMM  ♦  l 
ZERO  ■  0.0*0 

lh :  cors?;Lcxsi5 til 

XPTR  »  X*T 
XPTI 

CXPT  .  _ 

CS3  •  COEXP(CXPT) 

CWJM  ■  CS2*CS3 
CSA  »  C$l*CSl*CSl*CSl 
CS5  »  OCMOLXI AL .ZERO) 

CS6  «  CS5*CS5  *Cs5  *CS5 
COENOM  -  CS4  -  CS6 
CRESLT  ■  CNUM/COENOM 
U  -  OPEAL (CRESLT) 

V  »  01  MAGIC  RE  SLT) 

RETURN 
ENO 


DCMPLX  IXPTR,XPTI) 


AL  ■  1.2  50- 1 
T  -  8.D+0 
AA  -  0.3 
Xp  -  -4.99 
yF  -  2.3 
MMMMM  *  24 
LOSC  -  5 

SUM  ■  1.04169D+0 
SUMA  •  2. 18326D+0 
TIME  »  0.933140-1 
ANALYTICAL  ■  1.04169D+0 
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Case  62. 


£ 


- 1 


v<+«  + 


-  X(k  - 


SUBROUTINE  VALUE  IX,  V,  T,  AL ,  BE,  U,  V,  MMMMM,  KC  ALL) 

THIS  IS  "AELE  ET°Y  NUM9£0  62  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 

IMPLICIT  *"•  -  -  -  - 


,  __  IT  REAL*8  C A,8,D-H,3-ZJ 

I  Mol  IC  IT  r.3**PLr:X*l6  I  C  ) 

OI.MEK'IlCN  ME  ACL  (5  I 
NMHMM  ■  MMMMM  ♦  I 
I ERO  ■  0.0*0 

OCMPLX (X ,Y ) 

OC MPLX(AL,ZERO) 

DCMPLX (BE, ZERO) 

XPTR  «  X*T 
XPTI  -  Y*T 

CXOT  ■  OCMOLX  (XPTR.XPTII 
CNUM  •  COEXPICXPT) 

C$4  »  CS1  ♦  CS2 

cIll”cpOw£R  (!s4,CS6, 5.0-1 ,1  ,NCALU 
CALL  roOWER  <C$5,CS7,5.D-1,2,NCACL) 
COENOM  »  CS6  ♦  CS7 
CRESLT  ■  CNUM/COcNOM 
U  ■  OREAL(CRESLT) 

V  •  DIMAG(CRESLT) 

RETURN 

END 


AL  »  1.25D-1 
BE  -  -1 .250-1 
T  ■  8.EH-0 
AA  -  0.3 
Xp  -  -4.52 

yP  •  2.2 

MMMMM  -  23 
LOSC  »  5 

SUM  ■  1.17209D-1 
SUMA  -  6. 50880D-1 
TIME  »  0.13901D+0 
ANALYTICAL  «  1.17209D-1 
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Case  63 • 


£ 


- 1 


M  - 


Vi 


USRIUTIVE  V4LJE  (X,  Y,  T,  4L,  U,  V,  MMMMM,  NCALL) 

HIS  IS  ‘'‘ABLE  ENTRY  N'JMeER  63  OF  LAPLACE  TP  AN  SFORMS  BY  SPIEGEL 
-  Bjg-jH.O-2) 


M»LIC1T  b  EAL»8  (A, a 
‘  ldlT  COMPLEX* 16 


i  ppl 


DIMENSION  NCALH5) 

MHMHM  m  HMMMM  ♦  1 

ZERO  *  O.O+O 
CS1  -  DCNOLXI  X»Y) 

CS2  «  OCNPLXI AL.ZEROI 
XPTR  «  X** 

CXPT  M  OCMPLXIXPTR.XPTI ; 

CNUM  «  COEXPICXPT) 

C  S3  ■  CSl  ♦  CS2 

CALL  C POWER  CCS3, CS4, 5.0-1, l.NCALLI 
COENOM  »  CSl*CS4 
CRE5LT  •  CNUM/COEMON 
U  »  DREAL(CRESLT) 

V  »  DI MAGICRESLT) 

RET'JPM 

ENO 


AL  »  1.250-1 
T  ■  8.D+0 
AA  -  0.2 
Xp  ■  -4.64 

yF  ■  2.2 

MMMMM  »  23 

LOSC  -  5 
SUM  »  2. 383520*0 
SUMA  ■  2.647360*0 
TIME  •  0.852020-1 
ANALYTIC \t  «  2.383520+0 


155 


Case  6k. 


£  [/*<—)]  — 


«*  trtVZl 1 

Vi 


SUBROUTINE  VALUE  (X,  Y,  7,  AL,  U,  V,  MM**MM,  NCALL) 

THIS  IS  TABLE  ENIRY  NUMBER  64  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REAL*B  < A, 3, C-h,0-Z ) 

Implicit  complex* i6  <cj 


OIMENS  ICN  ‘!C.  ALL  (51 
MMMMM  •  MNMMM  «.  l 
ZERO  «  0.0*0 
CSX  *  OCMPLX(X.Y) 

CS2  »  DCMPLXIAL,Z£ROJ 
XO’-R  a  X*T 
XPTI  «  Y*T 

CXPT  »  DCMOLX(XPTR,XPTI) 

CNUM  ■  COEXP(CXPT) 

CALL  CPOWER  (CS1.CS3, 5.0-1, l.NCALU 
CSA  *  C  SI  -  CS2 
COENOM  >  CS3*CS4 
CRESLT  «  CNUM/ COE  NOM 
U  ■  OREAL(CP.ESLT) 

V  »  DIM  AG( C°  ESLT  J 

RETURN 

END 


AL  «  1.25D-1 
I  •  8.D+0 
AA  *  0.4 
Xp  -  -*-1* 

yp  *  2.2 

KMMMM  ■  23 
LOSC  -  5 

SUM  -  6.47907D+0 
SUMA  *  8.568020*0 
TIME  -  0.975520-1 
ANALYTICAL  ■  6.479070+0 
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Case  65. 


SUBROUTINE  VALUE  (X,  Y,  T,  AL,  BE,  U,  V.  MMMmh,  NCALL) 
THIS  IS  tABLE  ENTRY  NUMBER  65  3F  LAPLACE  TRAN»ORM$  BY 
IMPLICIT  REAL*8  (A,B,C-h,3-Z) 

Implicit  complex* ia  «c i 
01 MENS  I  ON  NCALL (5) 


♦  1 


:R0) 
ERO ) 


MMMMM  ■  MMMMM 
ZERO  *  0.3*0 
CSI  »  OCMPLXIX.YJ 
CS2  -  OCMPLXIAL.ZE 
CS3  »  OCM»LXI BEtZE 
xr*  ■  x*t 
XPTI  -  Y*T 

CXPT  »  OCMPLXI X® TR  *XPTI  ) 

CMJM  »  COEXPICXPT  I 
C$4  •  CSI  -  CS2 

CALL  CPOWEP  (CS4.CS5, 5.0-1, l.UCALL) 
CDENOM  -  CS3  +  CS 5 
CAESLT  »  CNUM/C3EN0M 
U  -  OREAL(CRESLT) 
v  -  oimagicreslt; 

*£TJRN 

EM) 


AL  -  1.25D-1 

HE  •  3.535533906D-1 

T  «  8.0+0 

AA  »  0.3 

Xf  -  -4.54 

yF  -  2.2 

MMMMM  ■  23 

LOSC  ■  5 

SUM  •  1.31286D-1 

SUMA  •  9.05849D-1 

TIME  •  9.01160-2 

ANALYTICAL  ■  1.31286D-1 
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SPIEGEL 


AL  *  1.25D-1 
T  •  8.D+0 
AA  «  0.2 
Xy  •  -4.64 

yF  -  2.2 

MMMMM  -  23 
Lose  «  5 

SUM  «  7.65198D-1 
SUMA  -  1.22916D+0 
TIME  ■  0. 91286D-1 
ANALYTICAL  =  7.65198D-1 
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Case  67 


_  A 


AL  «  1.25D-1 
T  -  8.D+0 
AA  «  0.3 
Xp  -  -4.99 
yF  -  2.3 
MMMMM  -  24 
LOSC  -  5 

SUM  ■  1.26607D+0 
SUMA  «  2.31076D+0 
TIME  -  0.10291D+0 
ANALYTICAL  ■  1.26607D+0 
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Case  68 


•■V.(at) 


SUBROUTINE  VALUE  (Xt  V.  T,  AL.  NPOWER «  U*  V.  MMMMM,  NCALL) 

THIS  IsTABLE  ENTRY  NliMR£&  68  OF  LAPLACE  TRANSFORMS  eY  SPIEGEL 
IMPLICIT  REAL  *8  (  A.3.0-HO-Z) 

IMPLICIT  C0mplEX*16  1C) 


IMPLICIT  REAL *8  < A.3,D-H,0-Z) 

IMPLICIT  COMPL  EX*  16  1C) 

DIMENSION  NCALLI5) 

MMNMN  •  MMNMN  4-  X 

ZERO  ■  O.DtO 
CSl  •  DCMPLXIX.Y) 

CS2  ■  OCMPLXl AL.ZERO) 

XPTR  -  X*T 

CXPT  «  CCM’LXIXPTR.XPTI ) 

CS3  «  <CSI*CS1J  t  ( CS2*CS2) 

CALL  CPOWER  IC$3,CDENO.5.0-ltl  .NCALL) 
CS*  »  (CDENOM  -  CS  L ) **N POWER 
CNUM  »  CS4*CDEXP( CX“T) 

CRESLT  ■  CNUM/CObNOM 
U  •  ORE ALICRESLT ) 

v  ■  oi magicreslt) 

RETJRN 

END 


AL  -  1.25D-1 
N POWER  -  3 
T  »  8.D+0 
AA  «  0.4 
Xp  «  -2.16 


MMMMM  *  17 
LOSC  «  1 

SUM  -  3.82096D-5 
SUMA  *  3.82096D-5 
TIME  *  0.81224D-1 
ANALYTICAL  •  3.82090D-1 
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<c«v>oooo 


Case  69 


-X  (.  -  y^nr^).  a>_j 

oC  l  v'72-“'  J 


AL  «  1.25D-1 
N POWER  =  3 
T  *  8.D+0 
AA  *  O.A 
Xp  -  -2.16 
yp  -  1.6 

MMMMM  =•  17 
LOSC  *  1 

SUM  =  U.32979D-5 
3UMA  *  U.32979D-5 
TIME  *  0.74672D-1 
ANALYTICAL  *  U.32977D-5 
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+  26) ) 


v;--_  - . 


Case  70. 


ItfPrS^SL^ENT&’NW*  ft’oFi^Ulk  TRANSFORMS  '’BY  SPIEGEL 
IMPLICIT  REM.  *8  ( 4,3 *D-H*0-Z ) 

IMPLICIT  C0‘1PLEX*16  (C) 

DIMENSION  NC4LL15) 

HMMMM  »  mmmmm  •*•  1 

ZERO  *  O.D+O 

CS1  -  DCMPLX(X.Y)  _ 

CS2  *  DCMPLXI 4L.ZESO ) 

CS3  *  DCMPLXiaE.ZERO) 

XPTR  «  X*T 
x  PT  I  »  Y  *T 

CXPT  «  DCMOLXIXPTRtXPTI) 

cciiL*ci.y»cf5u:c®^;!;i-L,i,«ALL. 

^num**  ' CD^X®( C S5 ) *CDE XP?CX°T ) 

CRESLT  ■  CNUM/ CDfcNOM 
U  «  DREAL(CRESLT) 

V  *  DIMAGI  CRESLT) 

RETURN 

END 


AL  =  1.25D-1 
BE  *  -1.25D-1 
T  =  8.D+0 
AA  =  0.2 
Xp  - 

yp  -  2.2 
MMMMM  *  23 
LOSC  -  5 

**  SUM  »  7.72O80D-1 
SUMA  »  1.23^85D+0 
TIME  -  0.82914D-1 
ANALYTICAL  ■  7.  72088D-1 
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Case  71 


<  <  6 


AL  -  1.25D-1 
BE  «  -1.253-1 
T  -  8.3+0 
AA  -  0.2 
Xp  «  -4.64 
yF  -  2.2 

MMMMM  *  23 

Lose  -  5 

SUM  *  7.65252D-1 
SUMA  -  1.246903+0 
TIME  ■  0.866583-1 
ANALYTICAL  =*  7.65252D-I 
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1 


Case  72. 


-1  _ i _  — 

JP  *  l  (^+7»pjj  —  a 


?2'o^'u»L»Tfl>AS5ftWs  BY  SPIBOEL 
IMPLICIT  REM.  *8  ( A,3,D-H,0-Z) 


IMPLICIT  C0MplEX*16  (C) 
DIMENSION  NCALL15) 

MMMMM  ■  M“MMM  ♦  1 


ZERO  »  O.O+O 
XPTR  «  X*T 
X  PTI  *  Y*T 

CXPT  *  OCMPLX<XPTR.XPTI) 

CNUM  »  COEXOICXPT) 

X  REAL  -  <X*X)  -  (Y*Y)  ♦  |AL*AL  ) 

XI  MAG  «  2.0+0*X*Y  .  „ , 

CS1  *  DCMOL  X(  XREAL.XIMAG) 

CS2  *  CSI*CU*CS1 _ ... 


CPQWER  (CS2 .  COE  NCR  *  5.0-1.  1.  NCALLI 


&RESLT  2  C‘IUm/Co£nOM 
U  «  OREAL(CRESLT) 

V  -  0  IMAGl  C  RE  SC  TI 

RETURN 

END 


AL  -  1.25D-1 


T  •  8.D+0 


AA  »  0.4 


Xy  -  -4,44 


yp  “2.2 


MMMMM  *  23 


LOSC  »  3 


SUM  *  2.81632D+1 


SUMA  «  2.«1715D+1 


TIME  ■  0.11476D+0 


ANALYTICAL  *  2.81632D+1 


^,OF'u<LKr-MyLn«is  »  SP.KEL 

1*  E' 


Implicit  ;omplex*16 

DIMENSION  NCALLI5) 

MPMMM  -  MMMMM  ♦  1 
ZERO  *  C.D  +  O 
CR1  «  OCMPLX(X.Y) 

X PTR  «  X*T 
XPTI  «  Y*T 

CXPT  ■  0CMPLX<X9TR,X»TI» 

C$2  -  COEXP  ICXPT) 

CNUM  "  CS1*CS2 

X REAL  »  <X*X)  -  (Y*Y)  ♦  (AL*AL) 
XJMAG  *  2.0+0*X*Y 
;S3  *  DCM°L  X (XRE4L*XIMAG) 


Se/cWitfW,  =0E«O.  5.0- 
CRESLT  *  CNUM/CDENOM 


1,1,  NCALL) 


0  »  DREAU  CRESLT) 
V  ■  3 1  MAGIC  RE SLT) 
RETURN 
END 


AL  =  1.253-1 
T  *  8.D+0 
AA  «  0.4 
Xp  *  -4.89 

yF  *  2.3 

KMMMM  *  24 

LOSC  «  4 

SUM  -  6.12158D+0 

SUMA  «  8.34998D+0 

TIME  -  0.12620D+0 

ANALYTICAL  «  6.12158D+0 
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Case  74. 


/,<•«)  - 


VALUE  (X«  Y«  T,  AL»  U,  V*  MMMMM,  NCALL) 

THIS  IS  TABLE  E1TRY  NUMBER  74  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  R  E  AL  *  9  (  At  9  •  D-Ff  O-Z  ) 

IMPLICIT  C0moLEX*16  (Cl 
OIMENS  ION  '1C»LL(  5) 

NMMMM  »  MMMMM  ♦  1 

ZERO  ■  O.O+O 
CSI  -  OCMt>LX(X,YI 
XPTR  «  X*T 
XPTI  *  Y*T 

CXPT  »  OC«PLX(XPTRt XPTI ) 

CS2  *  COEXPICXPT) 

C  S3  *  CS1*CS1 
CNJM  *  CS2*CS3 

XREAL  -  ( X*X)  -  (Y*Y)  ♦  ( AL*AL ) 

XIMAG  -  2.D+0*X*Y 
CS4  *  OCMPLXIXREAL.X IMA3I 
C  S5  ■  C  S4*CS4*CS4 
CALL  CPOWER  (CS5.  COENOM, 

CRESLr  -  CNUM/COE  N0M 
U  ■  OREAL(CRESLT) 

V  -  DIM AG(CRSSLT) 

RETURN 
ENO 


5.0-1,  1,  NCALL) 


All  -  1.25D-1 
T  *  8.D+0 
AA  »0.4 
Xj.  ■  -4.89 
yp  -  2.3 

MMMMM  ■  24 

LOSC  «  5 

SUM  -  3 .251470-1 

SUMA  «  3.725U7D+0 

TIME  ■  0. 10772D+0 

ANALYTICAL  «  3-251470-1 
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mvn  'n.rt'-Fwtfsife.fe  Aarf^sM 

iw»li:it  RE4L*a  <  ft. a.  o-h,o-z) 

Implicit  C1M°tEX*16  <  C) 

01  HENS  I  ON  NC  ALL (5  > 

MMMMN  •  MMNMM  ♦  1 
ZERO  »  0.0+0 
XPTR  -  X*" 

xpti  •  y*t 

CX®T  -  OC.MPLXIXPTR,  XPTII 
C  MJM  »  COEXPICXPT) 

XP6AL  •  (X*XJ  -  I  Y*  Y I  -  !AL*4L) 

XI MAG  •  2.0+O*X*Y 

cjl  .  0CMOL XI XREAL.XIMAG) 

CALL*CPOWER5 (CS2t  COENOMt  5.0-l»  1*  NCALLI 
CRESLT  -  CNUN/COENOM 
•J  •  O'*  EAL I  CRESLT) 


fRMS  BY  SPIEGEL 


NCALLI 


01NAGICRESLT I 
URN 


AL  -  1.25D-1 
T  ■  8.D+0 
AA  ■  0.5 
Xp  -  -U.34 

yF  -  2-2 

KKMMM  *  23 
LOSC  ■  3 

SUM  •  3.61702D+1 
SUMA  ■  3.638U6D+1 
TIME  -  0. 110090+0 
ANALYTICAL  «  3.6X7020+1 
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AL  =  1.25D-1 
T  »  8.D+0 
M  *  0.5 

Xp  «  -4.79 

yp  -  2.3 

MMMMM  *  24 
LOSC  -  4 

SUM  »  1.01295D+1 
3UMA  -  1.52218D+1 
TIME  ■  0.13250D+0 
ANALYTICAL  -  1.012S5D+1 
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Case  77- 


;,(•«)  +  •»/,(«<) 


C 


BY  SPIEMl 

IMPLICIT  REAL*8  l O-Hb^-Z ) 

IMPUCIT  C0>1PLEX*I6  (Cl 
DIMENSION  NC  ALL ( 5) 

MMMMM  ■  ♦  l 

ZERO  *  0.0*0 
CSl  *  OCM°L  X(  X.  Y) 
xrR  *  x*t 

cxpl  «  OcJplXCXPTP.,XPTI) 


•XPT  «  OC*PLX(XPTP ,XPTi: 

m :  cs!2cs? 


XR6AL"»CU*Xi3-  (Y*Y)  -  (AL*AL) 
i|?S  DCMPLX*(XR  EAL»  X  IN AG  I 

CALL*CPOWEP Stcsst  CDENQM,  5.0-1.  1*  NCALL) 
CRESLT  —■ 


,1  1  V.KUWC''  (  W  J 

t^mr 


RETURN 
END 


AL  »  1.25D-1 
T  =  8.D+0 
AA  -  0.5 
Xy.  «  -4.79 

yp  -  2.3 

MMMKM  *  24 
LOSC  -  5 

SUM  -  1.33122D+0 
SUMA  *  7.71705D+0 
TIME  =  0.1I807D+0 
ANALYTICAL  *  1.83122D+0 


169 


Case  78 


CO.tySI 

vr, 


SUBROUT IME  VALUE  <X,  Y,  T,  al.  U,  V,  MMMMM.  NC ALL) 

THIS  IS  TABLE  ENTRY  NUMBER  81  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  PEAL*fl  (A.P.O-H ,0-ZJ 
IMPLICIT  C0MPLEX*16  <  C I 
DIMENSION  NCALLI5) 

MMMMM  m  MMMMM  1 

ZERO  ■  O.D+O 
CSI  -  OCMPLXJX,  Y) 

CS2  »  DCM9LXIAL.ZER0) 

CS3  •  -CS2/CSI 
XPTR  ■  X*T 
XPT I  -  Y*T 

CXPT  »  0CM°LX (X  PT  R,X  °T I ) 

CNUM  ■  COEXP(CSA)  *COEXPCCXPT) 

CALL  CPCWER  CCS l, CCcNOM, 5.0-l» l.NCALL) 

CRE$LT  *  C  NUM/COENOM 
U  »  DREALCCPESLT) 

V  »  DIMAGCCRESLT) 

R  ETUR  N 
END 


AL  «  1.25D-1 
T  =  8.D+0 
AA  »  0.2 
Xy  *  -4.64 
yp  -  2.2 

MMMMM  *  23 
LOSC  »  5 

SUM  -  -8.300930-2 
SUMA  ■  4.47194D-1 
TIME  ■  0.96148D-1 
ANALYTICAL  *  -8.30093D-2 
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Case  79 


Vj  -1  — 

£  - 


-  tin  ZJTt 

—  vr* 


SUBROUTINE  VALUE  (X,  Y,  T,  AL,  U,  V,  MMMMM,  NCALU 

THIS  IS  TABLE  ENTRY  NUMBER  82  OF  LAPLACE  TRANSFORMS  BY  SPIESEL 

IMPLICIT  RE AL*8  ( 4. 3 , C-htO-2 ) 

Implicit  co«plex*16  <c> 

DIMENSION  NC ALL (51 
MMMMM  m  MMMMM  ♦  1 
ZERO  -  O.O+O 
CSl  ■  DCMPLXIX.Y) 

C  $2  »  DCMOLX(AL.ZERO) 

CS3  «  -CS2/CSI 
XPTR  -  X*T 
XPTI  ■  Y*T 

CXPT  »  CCMPLX(XPrR,X«>TI) 

C  MJM  -  COE X PICS 3) *CDEXPICXPT ) 

CS4  ■  CS1*CS1*CS1 

CALL  C POWER  (CS4, CCENOM.S.D-l,  l.NCALLI 
CRESLT  *  CNUM/COENOM 
U  «  DREALICESLTI 
V  ■  DIMAG(CRESLT) 

RETURN 

END 


AL  -  1.2 50-1 
T  «  8.D+0 
AA  >  0.2 
Xy  ■  -4.64 

yF  •  2.2 

MMMMM  -  23 

LOSC  -  5 

SUM  -  1.45103D+0 

SUMA  ■  1.82676D+0 

TIME  -  0. 10018D+0 

ANALYTICAL  *  1.45103IH0 
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Case  80 


•  Zi,u  «•  oJ'SRicli’T^rsSSi  »*M>  Ieoel 

IMPLICIT  REM *8  (  1,3.0-H,0-Z) 

IMPLICIT  COMPLEX* 16  <C) 

0  IMENSI ON  NC  ALL (5) 

MMMMM  ■  ♦  1 

ZERO  »  O.D+O 
CS1  -  g^MPLXIX,  Y) 


CS2 

XPTR  « 
XPTI  « 
,XPT  » 
,NUM  « 
COENI* 
CRESLT 


MPLX (AL. ZERO) 


A 

s 


X*T 

r*T 


DCM»LX(XOT  R,XP-n. 
CDEXP(CS3)*C0E) 


. . . . XPICXP 

CS1**I  MOOWEO  *  1) 

_  CMUM/COENOM 

U  ■  OPE ALI CRESLT) 

V  ■  0IM»5(CRESLT) 

RETURN 

END 


AL  -  1.25D-1 

NPOWER  ■  3 

T  ■  8.D+0 

AA  -  0.3 

Xp  -  -U.il 

yp  ■  2*1 

MMMMM  -  22 

LOSC  «  U 

SUM  -  6.60189D+1 

SUMA  ■  6.60189D+1 

TIME  »  6.9758D-1 

ANALYTICAL  -  6.60139D+1 
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SUBROUT IME^V  ALU CtU.^Y.^Tj  'L ,  ^U,  ^j^mmmm^NCALL) $  ^  $p  j  EGEL 
IMPLICIT  0  EtL *8  ( A.B , D-rt,3-Z ) 


I(*PLICIT  COMPLEX*  16  (C) 

0 IME'l  SI  ON  NCALL<  5) 

MMMMM  ■  MMMMM  ♦  1 

ZERO  ■  0 .D+O 
CS1  »  OCM°LX(  Xf  Y) 

CS2  -  OCMP  LX ( AL • ZERO ) 

XPTR  *  X*T 
XPTI  *  Y*T 

CXPT  *  OCMPLX  (X°TR,XPT  I ) 

CALL  C POWER  (CS1 ,COENCM, 5.0-1 .1 .MC4LLI 

CS3  *  -CS2*CDEN0M 

CNUM  ■  CDEX  P(  CS3  )  *CQEXP  t  CXPT  ) 

CRESLT  *  CMUM/COENOM 
U  -  DR EAL( CRESLT) 

V  »  01 M AG t CRESLT ) 

R6TURM 

EMO 


AL  *  1.25D-1 
T  «  8.D+0 
AA  »  0.2 
Xp  ■  -4. 61+ 

yF  -  2.2 
MMMMM  *  23 
LOSC  -  5 

SUM  -  1.99374D-1 
SUMA  •  6.27444D-1 
TIME  *  0.93470D-1 
ANALYTICAL  *  1.99374D-1 
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Case  82. 


ly/7t* 


TH?S°i*S  *'*'A8LE*"2NTRY f  NUM3E3  85*  PF*  LA^LACE^T&ANSFORMS  BY  SPIEGEL 
ILLICIT  5.E1L*8  (  A« 8  * D-H ,3-Z ) 

IMPLICIT  C-l^PLEX*!^  ICJ 
DIMENSION  NC ALL ( 5 ) 

MMMMM  »  MMMMM  L 
ZERO  »  0.0*0 
C£I  •  OCMPl_X{  X»YI 
CS2  ■  OCM°LXt 4L. ZERO) 

XPTR  *  X*Z 
XPTI  ■  Y*T 

CXPT  »  DCMPLXIXPTR.XPTI)  ,  -  - 

CALL  CPOWER  l CS 1 »  Ca  3  »  5 • C-l t  X  *  NC  AL  L ) 

CPESLT  *SCDEXP(CS4J*C0EXP( CXPT J 
U  -  OREAUCRESLTI 
V  ■  DIM AG(CRESLT) 

RETURN 

EM) 


AL  -  1.25D-1 
T  *  8.D+0 
AA  -  0.2 
Xp  -  -5.09 
yF  -  2.3 
MMMMM  *  24 
LOSC  ■  6 

SUM  «  1.55761D-3 
SUMA  -  3.03773D-1 
TIME  -  0.10197D+0 
ANALYTICAL  ■  1.55761D-3 
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Case  83. 


c 


f  HIS?  £S  ^PLE^ENT  AY  lN£^n_f  ^  ^  ’  L  APL  ACEMtJ  AN  SFORM  S  BY  SPIEGEL 


IMPLICIT  REAL*8  ( At B t 0-H» Q-Z I 
IMPLICIT  COMPLEX*  16  I  C) 

01 MENS  ION  MC  ALL ( 5  ) 

MMMMM  «  MMMMM  ♦  I 
ZEPO  *  0.0*0 
OKE  -  l.D*0 
CS1  *  OCMPLXIX.YI 
CS2  *  DCMOLXI AL tZEPTl 
C$3  «  OCMPLXIONE,  ZERO) 

XPTR  «  X*T 

11 

V  «  DIMAGICRESLT) 

RETURN 

END 


AL  ■  1.25D-1 

T  -  8.D+0 

AA  *  0.2 

Xj.  -  -4.64 

yF  «  2.2 

MMMMM  »  23 

LOSC  ■  5 

SUM  -  2.49298D-2 

SUMA  -  8.07480D-2 

TIME  -  0.93158D-1 

ANALYTICAL  -  2.49292D-2 
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Case  84.  — » 


trie (a/Zy/t  ) 


SUBROUTINE  VALUE  (  X,  Y,  T,  AL,  U,  V,  MMMMM,  NCALL) 

THIS  IS  "ABLE  ENTRY  NUMBER  87  "*F  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLI  Cl T  R EAL*B  ( A ,3 ,D-H,0-Z) 

IMPLICIT  COMPLEX* lt>  (C) 


IMPLICIT  COMPLEX*l«> 
OIMENSION  NCALLI  5  I 
MMMMM  «  m^MMM  >  1 

ZERO  =  O.O+O 
CSl  -  DCMPLXU.Y) 

C  S2  *  OCM°LX( AL.ZSRO) 


X  PTR 
XPTI 
CXPT 
CALL 
CNUM 
COENOM 


X*’ 

«  Y*T 

*  OCMPLX(XPTR,XPTI) 

CPONER  (CS1, CS4, 5.0-1, 1, NCALL) 

CDEXP(-CS2*CS4)*C0EXP(CXPT  ) 

_  •  r$l 

CRESLT  *  CVUM/ COENOM 
U  -  OREAL(CRE  SLTI 
V  x  DIM AGJCRESLT ) 

RETURN 

END 


AL  =  1.25D-1 
T  m  8.D+0 
AA  *  0.2 
X?  »  -4.64 

yp  *  2*2 

mmmmm  *  23 
LOSC  ■  5 

SUM  «  9.75070D-1 
SUMA  *  1.35850D+0 
TIME  •  0.99658D-1 
ANALYTICAL  =  9.75071D-1 
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SUBROUTINE  VALUE  IX,  Y,  t ,  AL ,  RE,  U,  V,  MMMMM,  NC  ALL ) 

THIS  IS  TABLE  ENTRY  NUMBER  88  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 


IMPLICIT  REAL *8  (  A,  8,  C-F..O-Z  J 
IMPLICIT  COMPLEXM6  (CJ 
DIMENSION  N CALL ( 5  J 
MMMMM  *  MMMMM  +  1 

ZERO  «  O.O+O 
CSX  »  OCMPLXIX,  Y) 

CS2  ■  OCMPLX (  AL,  ZERO! 

XPTR  «  X*T 
XPTI  *  Y*T 

CXPT  *  OCMt»LX(X°TR,X?Tl) 

C  S3  ■  OCMPLXIRE.ZER  T) 

CALL  CPOWER  ICS  I, CS4,  5.0-1,  l.NCALLJ 

CNUM  «  COEX P (-CS2+C3 4 ) *COEXP ( CX PT  ) 

CS5  *  C  S3  ♦  CS4 

COENCIM  ■  CS4+CS5 

CRESLT  «  CNUM/COENOM 

U  «  DREALICRESLTJ 

V  *  DIMAGICPESLTI 

RETURN 

END 


AL  «  1.25D-1 

BE  »  3.457408906D-1 

T  ■  8.D+0 

AA  ■  0.2 

Xy  •  -4.64 

y,  *  2  2 

MMMMM  «  23 
LOSC  -  5 

SUM  «  4.27375D-1 
SUMA  -  7.5575D-1 
TIKE  -  9. 92 940-2 
ANALYTICAL  ■  4.27375D-I 
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Case  86 


« 


c 


SUBROUTINE  VALUE  (  X  »Y  , T ,AL  »8E  ,U 
ThIS  IS  TABLE  ENTRY  NUMBER  90 
IMPLICIT  REAL*8  ( A,9,0-H,T-Z) 
IMPLICIT  COMPLEX*  lo  (C) 
DIMENSION  NCALLI5I,  MCALL ( S I 
MMMMM  «  MMMMM  ♦  1 
ZERO  =  0.0*0 
CS1  ■  DCMPLXIX.Y) 

CS2  «  OCMPLXIAL.ZERO) 

C S3  «  DCMPLXI 6E.ZER0) 

XPTR  «  X*T 
XPTI  «  Y*T 

CXPT  ■  DCMPLXI  XPTR,  XPTI) 

CS4  «  CSI  ♦  CS2 
CS5  «  CSI  ♦  C S3 
CS6  ■  CSA/CS5 


DREAL  fC RESLT) 
V  -  DIMAG(CRESLT) 
RETURN 
END 


MCALL) 


) 

BY  SPIEGEL 


AL  -  1.25D-1 
BE  -  2.5D-1 
T  *  8.D+0 
AA  *  0.2 
Xr  -  -D.21 

yP  •  2-1 

MMMMM  *  22 
LOSC  ■  5 

SUM  •  -2.906800-2 
SUMA  ■  1.017«OD-1 
TIME  •  7.00960-2 
ANALYTICAL  *  -2.906800-2 
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Case  87. 


"1  I  la  [<«*  +  «•)/•*]  — 

L^J  - 


SU8>"1UTI'IE  VAL'Jc  (Xf  Y,  Tt  4L  *  Ut  Vi  MMMMM,  MC4LL) 

^  LAPLAC£  7RA^P0RMS  BY  SPIEGEL 

J  PPL  I C I T  CO‘'PL  EX*  16  I  C  J 


IMPLICIT  COMPLEX *16  1 C)  ’  '  “ 

DIMENSION  MCALL  (5) 

MMMMM  ■  MMMMM  ♦  1 
ZERO  ■  O.D+O 
CSl  -  OCMPL X( Xf  YJ 
CS2  »  OCMPLX< AL.ZERO) 

CS3  -  HCS1*CS1)  ♦  (  CS2*CS2) )/(  CS  2*CS2) 
XPTR  ■  X*T 
XPTI  ■  Y*T 

cxpt  ■  oc«plx ixpri.xPT n 

CALL  CLOG  (CS3,  C  S4 ,  1,  MCALLJ 
CNUM  -  CS4*CDEX?<  CXPT ) 

£ S5  «  DC mplx <2. O+O, ZERO) 

CRESLT  «  CNUM/CCS5+CS1) 

U  •  ORE  AL ( CRESLT  J 
V  «  OIMAGICRESLTJ 
RETURN 
END 


AL  ■  1.25D-1 
T  ■  8.D+0 
AA  «  0.2 

Xy  •  -5.09 
yF  •  2.3 

MMMMM  -  24 
LOSC  -  5 

SUM  »  -3*37404D-1 
SUMA  ■  1.27603D+0 
TIKE  •  8.8244D-2 
ANALYTICAL  ■  -3.37404D-1 


Case  88 


!'  '"^"’"5^’  l"J  ■'WWIW 


El(«0 


SUBROUTINE  VALUE  (X,  Y,  T,  ALt  U,  V,  MMMMM,  MCALL) 

THIS  IS  TABLE  ENTRY  NUMBER  92  OF  LAPLACE  TRANSFORMS  BY  SPIEoEL 
IMPLICIT  S  EAL»8  ( A,B, O-H.O-Z) 

IPPUCIT  C0MPLEX*16  (Cl 
DIMENSION  CALL  (5) 

MMMMM  *  MMMMM  ♦  1 
ZERO  »  G.Ot-O 
CS1  ■  OCMPLX(X.Y) 

CS2  •  DCMPLX(AL,ZER0J 
C  S3  ■  (CS1  ♦  CS2J/CS2 
XPTR  «  X*T 

•  DCMPLX(XPTR,XPTI» 

CALL  CLOG  ( CS  3*  CS4,  1.  ‘CALL ) 

CNUM  •  CS4*CD£XP(  CX“T> 

CRESLT  ■  CNUM/CSI 
U  «  OPEALICRESLTI 
V  -  DIM AG(CRESLT) 

RETURN 

END 


AL  ■  1.25D-1 
T  »  8.D+0 
AA  -  0.2 
Xp  -  -5.09 
yF  •  2.3 

MMMMM  ■  24 
LOSC  •  5 

SUM  •  2.19384D-1 
SUMA  •  1. 55094 OtO 
TIME  -  7.2124D-2 
ANALYTICAL  •  2.19384D-1 


Case  89 

Z'1 


-(!■»  +  r> 

y  a  Bulfr'l  rontUnt  » 


SUBROUTINE  VALUE  (X,  V,  T,  AL .  U,  V,  MMMMM.  MCALLJ 
THIS  IS  TABLE  ENTRY  NUMBER  96  OF  LAPLACE  TRANSFORMS 
IMPLICIT  SE'.L*8  <  A.B.O-H.D-Z) 

IMPLICIT  C0MPLEX*16  I C i 
DIMENSION  “CALL  (51 

MMMMM  m  MVMMM  +  1 

ZERO  ■  0.0*0 
CSi  »  DC MPL  X(  Xi  Y) 

.ALL  CLOG  (CSI 

>  X*T 

>  Y*T 

>  QCMOLXIXPTR.XPTn 

>  C52*C0cXPf  CXPT) 

CRESLT  -  CNUM/CSl 
U  ■  DREAL(COESLT) 

V  •  01  MAG(CRESLT) 

RETURN 
END 


II*  CS2t  1.  MC ALL ) 


AL  -  1.25D-1 
T  -  8.00 
AA  •  0.2 
Xp  -  -5.09 

yF  •  2-3 

MMMMM  ■  2d 
LOSC  -  5 

SUM  ■  -2.6566600 
SUMA  -  2.8060600 
TIME  ■  7.78180-2 
ANALYTICAL  •  -2. 6 5666  00 


BY  SPIEGEL 
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Case  90. 


u»i  +  r>*  -  i«* 

1  =  KulH'inudnt  =  .57TJ156. . . 


SUBROUTINE  VALUE  (X.  V,  T,  4 L *  U,  V,  MMMMM,  MCALU  ,  „„ 

THIS  IS  TABLE  ENTRY  NUMBER  97  OF  LAPLACE  TRANSFORMS  8Y  SPIEGEL 
IMPLICIT  REAL*8  (  A ,8 ^O-H , 0-Z) 


IMPLICIT  COMPLEX* 16 
DIMENSION  MC ALL  (5) 
mmmmm  »  MMMMM  *  l 
ZERO  «  0.0*0 
CSI  »  OC MPLXIX* Y) 

XPTR  ■  X*T 
XPTI  •  Y*T 

CXPT  »  DCMPLX  (XPTR,  XPTI) 

CALL  CL PG  (CSI,  CS2.  I,  m:4LL) 
CHUM  «  CS2  *CS2  *CDEXP ( CX°T ) 
CRESLT  -  CNUM/CSl 
U  »  OREAL(CPESLT) 

V  «  OIMAG(CRESLT) 

RETURN 

ENO 


AL  •  1.25D-1 

t  rn  8.0*0 

AA  •  0.2 
Xp  -  -5.09 
yp  -  2.3 

MMMMM  *  24 

LOSC  ■  4 

SUM  -  5.412890*0 

SUMA  »  5.434470*0 

TIME  -  7.73760-2 

ANALYTICAL  ■  5.412890*0 
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Case  91 


-  f  Y.  M 


IMPL  I C IT  cEM-*8  (  A,8»0-H,'W) 
iwucr  COMPLEX*  16  <Ci 
DIMENSION  NCALU  5),  MCALK5) 

MMHMM  ■  mm^mm  ♦  1 

ii?°.mOCMPLX(X,Y) 

DCM» LX ( AL i Z ERO I 

.  *  X*T 

XPTI  *  Y*T 

CXPT  «  0CMOLX (XPTR,XPTI) 

fa’cMy»*?ai,so4ss53ii.i.'.c»tLi 
fft-eispres*:  hv/fiw. 

lSi\-r  SS?«8%S^V 1 

U  «  3PE  AL(CAESLT| 

V  •  DIM AG<  CRESLT J 

RETURN 

ENO 


AL  •  1.25D-1 
T  •  8.D+0 
AA  -  0.3 
Xp  »  -4.99 
yF  «  2-3 

MMMMM  -  24 
LOSC  ■  5 

SUM  »  -1.386340-1 
SUMA  •  3.92274D+0 
TIME  ■  1.2126D+0 
ANALYTICAL  ■  -1.38634D-1 


183 


Case  92 


-ikc-r^]]  - 

[i*  *  » 


f  H.(«0 


This  is  known  as  the  Struve  function.  See  ref.  (1), 
page  496. 


SUBROUTINE  VALUE  (  X,  Y  ,  T ,AL  ,3  :  ,U  , V  .MMMMM.NC  ALL  ,MC ALU 
▼HZ S  IS  TAPL6  S “TRY  <12  *>*GE  2 64  FORESTS  AND  KAUFMAN 

JMPLICIT  5  EAL*d  (  4,3  ,D-H,0-Z) 

NPLICJT  COMPLEX* 16  t  C  J 
I MEN? I  ON  NCALLI5J,  MCALLI5) 

MMMMM  *  M**MMM  +  I 
ZERO  «  O.OfO 
C|l  »  DCMPLX(X.Y) 

C $2  ■  OCMP  LX( Ac  .ZERO) 


OIMEN? I  ON  MCALL1 5) .  MCALLC5) 

MMMMM  »  MMMMM  +  l 
ZERO  »  O.O+O 
C|i  «  OCMPLX(X.Y) 

C  $2  ■  OCMP  LX(  At., ZERO) 

XPTR  »  X*' 

XPTI  ■  Y*“ 

CXPT  •  OC"PLX<X?TR,XPTI) 

CS3  »  IIC?l*CSl)  ♦  (CS2*CS21 ) 

CALL  CPOWEF.  <CS3, COENO*. 5.0-1 ,1 .NCALL) 

CS4  •  ICOENOM  >  /c$l 

CALL  CLpG  (  CT4,  C?  a.  l  .  MCALL  ) 

ENUN  •  CS5*CDEXP(CXPT) 

RESLT  »  CriUM/cOEMOM 

uv :  htiWJUt’ji 

A  ET  JR  N 
ENO 


AL  -  1.25D-1 
T  »  8.D+0 
AA  -  0.3 
»  -4.54 
yr  »  2.2 
MMMMM  ■  23 
LOSC  •  4 
SUM  •  6.932440-1 
SUMA  •  9.27960D-1 
TIME  •  1.134 ID- 1 
ANALYTICAL  -  8.932440-1 
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Case  94 


SUBROUTINE  VALUE  (X , V ,T f AL.U, V . VMMMM ) 

THIS  TRANSFORM  WAS  TAKEN  FROM  R06ERTS  ANO  KAUFMAN 
IMPLICIT  PEAL *8  «A-H,  0-2) 

MPMMM  ■  MMMMM  «.  I 

TEM  1  ■  -X. 0*0*AL*AL*  X 

TEM2  ■  AL*AL*Y 

TEM3  «  4.0+  C*(  I  X*  X)  ♦  ( Y *Y ) ) 

TEM4  «  TEM  X/TEM  3 
~EM2 /T£M3 

|XP ( TEM4) *0C  OS  I  T£  MR ) 

_  §X°(TEM4)*0SIN(TEM5) 

_  ITEM6*X)  *■  (TE't7*Y) 

TEM9  -  (TEM7*X)  -  (TEM6* Y) 

TEM13  -  (X*X)  ♦  ( Y*Y ) 

TEM11  -  TgMg/xEMXO 


T|m.12 

TEM 

V  ■  iTEM14*T§MU) 
RETURN 
ENO 


TE“9/TE*U0 
M13  ■  OEXPIX«T)*OCOS(Y*T) 

14  -  OEXPl X*T) *0  SI N( Y*T ) 

( TEM X3*7tM XI)  -  (  TEM14*T=MX2) 
ITEM14*temi1,  ♦  (TEMX3*tEmi2) 


AL  ■  2. 50-1 
T  -  1.6041 
AA  »  0.2 
Xp  •  -3.04 
yT  •  1.8 
MMMMM  -  19 
LOSC  -  8 

SUM  *  7.65198D-1 
SUMA  •  3.98573D+0 
TIME  ■  6.3648D-2 
ANALYTICAL  ■  7.651980-1 


186 


Case  95 


* 


-at 

1  -  4at« 


SUBROUTINE  VALUE  I X  ,Y  ,T , AL.U.  V , 

THIS  TRANSFORM  WAS  TAKEN  FROM  ROBERTS  ANO  KAUFMAN 
IMPLICIT  R6AL*3  l A, 3 , 0-H, 0-Z) 

IMPLIC I T  COMPLEX* 16  1C) 

MMMMM  «  HMMHM  ♦  1 

ZERO  *  0. D>0 
CSl  ■  OCMPLXI  X,  Y) 

C$2  ■  OCMPLX (AL, ZERO! 

CS3  »  CSl  -  CS2 
CS4  ■  CSl  ♦  CS2 
XPTR  ■  X*T 
XPTI  «  Y*T 

CXPT  ■  DCMPLXIXPTR,XPTI ) 

CNUM  -  CS3*CS 3*CJEXP(  CXPT ) 

COENOM  *  C$l*CS4*CS4 
CRESL"  »  CNUM/COENGM 
U  -  ORE AL(CRESLT) 

V  *  DIMAG(CAESLT) 

RETJRN 

END 


AL  -  1.25D-1 
T  »  8.D+0 
AA  ■  0.2 
Xp  ■  -4.64 

yF  •  2.2 

MMMMM  ■  23 
LOSC  ■  5 

SUM  -  -4.715180-1 
SUMA  •  5.06864D-1 
TIME  -  6.9186D-2 
ANALYTICAL  -  -4.715J.8D-1 
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